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allocating  new  sample  points  in  regions  that  are  multi-modal  and  have  sensitive  response  versus  allocating 
sample  points  in  regions  that  are  sparsely  sampled.  In  the  second  research  thrust,  a  new  approach  for 
multiresponse  adaptive  DOE  is  developed  (i.e.,  MSFCVT).  Here  the  approach  from  the  first  thrust  is 
extended  with  the  notion  that  the  tradeoff  should  also  consider  all  responses.  SFCVT  is  compared  with 
three  other  methods  from  the  literature  (i.e.  maximum  entropy  design,  maximin  scaled  distance,  and 
accumulative  error).  It  was  found  that  the  SFCVT  method  leads  to  better  performing  metamodels  for 
majority  of  the  test  problems.  The  MSFCVT  method  is  also  compared  with  two  adaptive  DOE  methods 
from  the  literature  and  is  shown  to  yield  better  metamodels,  resulting  in  fewer  function  calls.  In  the  third 
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level  performance  evaluation  and  design  optimization  of  air-cooled  heat  exchangers. 
The  dissertation  is  divided  into  three  research  thrusts.  The  first  thrust  presents  a  new 
adaptive  DOE  method  for  single  response  deterministic  computer  simulations,  also 
called  SFCVT.  For  SFCVT,  the  problem  of  adaptive  DOE  is  posed  as  a  bi-objective 
optimization  problem.  The  two  objectives  in  this  problem,  i.e.,  a  cross  validation  error 
criterion  and  a  space-filling  criterion,  are  chosen  based  on  the  notion  that  the  DOE 
method  has  to  make  a  tradeoff  between  allocating  new  sample  points  in  regions  that 
are  multi-modal  and  have  sensitive  response  versus  allocating  sample  points  in 
regions  that  are  sparsely  sampled.  In  the  second  research  thrust,  a  new  approach  for 
multiresponse  adaptive  DOE  is  developed  (i.e.,  MSFCVT).  Here  the  approach  from 
the  first  thrust  is  extended  with  the  notion  that  the  tradeoff  should  also  consider  all 
responses.  SFCVT  is  compared  with  three  other  methods  from  the  literature  (i.e., 
maximum  entropy  design,  maximin  scaled  distance,  and  accumulative  error).  It  was 
found  that  the  SFCVT  method  leads  to  better  performing  metamodels  for  majority  of 
the  test  problems.  The  MSFCVT  method  is  also  compared  with  two  adaptive  DOE 
methods  from  the  literature  and  is  shown  to  yield  better  metamodels,  resulting  in 
fewer  function  calls. 

In  the  third  research  thrust,  an  approximation-based  framework  is  developed  for 
the  performance  evaluation  and  design  optimization  of  novel  heat  exchangers.  There 
are  two  parts  to  this  research  thrust.  First,  is  a  new  multi-level  perfonnance  evaluation 
method  for  air-cooled  heat  exchangers  in  which  conventional  3D  Computational 
Fluid  Dynamics  (CFD)  simulation  is  replaced  with  a  2D  CFD  simulation  coupled 
with  an  s-NTU  based  heat  exchanger  model.  In  the  second  part,  the  methods 


developed  in  research  thrusts  1  and  2  are  used  for  design  optimization  of  heat 
exchangers.  The  optimal  solutions  from  the  methods  in  this  thrust  have  44%  less 
volume  and  utilize  61%  less  material  when  compared  to  the  current  state  of  the  art 
microchannel  heat  exchangers.  Compared  to  3D  CFD,  the  overall  computational 
savings  is  greater  than  95%. 
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Chapter  1  Introduction 


In  this  chapter  we  present  the  motivation  behind  this  dissertation,  the  research 
objectives,  and  the  underlying  assumptions,  followed  by  a  brief  description  of  the 
research  thrusts. 

1.1  Motivation 

Real  world  engineering  systems  are  highly  complex  in  nature  and  often  involve 
multiple  disciplines  such  as  heat  transfer  analysis,  fluid  flow  analysis,  structural 
analysis,  control  systems  etc.  The  design  of  these  systems  is  more  complex  even 
considering  the  recent  advances  made  in  the  field  of  distributed  and  high  performance 
computing.  The  design  and  analysis  of  the  individual  disciplines  is  generally  carried 
out  with  the  help  of  complex  simulation  models.  An  example  of  such  a  complex 
simulation  is  a  Computational  Fluid  Dynamics  (CFD)  simulation  for  a  heat 
exchanger.  This  simulation  can  take  anywhere  between  a  few  minutes  to  several  days 
to  execute  and  converge,  not  to  mention  that  this  is  the  time  required  for  a  single 
simulation.  During  optimization  in  a  single  or  multidisciplinary  environment,  such  a 
simulation  may  be  invoked  by  the  optimizer  several  hundred,  even  several  thousand 
times.  Clearly,  carrying  out  a  design  optimization  with  such  a  complex  simulation 
tool  is  not  feasible.  Figure  1.1a  shows  this  conventional  optimization  procedure.  A 
complex  engineering  simulation  such  as  CFD  (Versteeg  and  Malalasekera,  1995; 
Anderson,  1995)  is  used  to  represent/model  the  performance  of  a  heat  exchanger  and 
is  coupled  with  an  optimization  or  other  analysis  tools  to  optimize  designs.  In  order  to 
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alleviate  the  computational  burden  associated  with  the  optimization  of  complex 
simulations,  researchers  have  proposed  a  technique  called  approximation-assisted 
optimization. 

Approximation  can  generally  be  summarized  as  a  two-phase  process:  (a)  design 
of  experiments  (DOE)  phase  and  (b)  development  of  a  metamodel.  In  the  DOE  phase, 
the  first  step,  the  complex  simulation  tool  is  run  for  a  set  of  systematically  chosen 
input  design  points.  In  the  second  phase,  the  inputs  design  points  along  with  the 
corresponding  simulation  results  (i.e.,  response  values)  are  used  to  develop  a 
metamodel.  A  metamodel,  which  is  also  referred  to  as  a  surrogate  model,  is  a 
simplified  numerical  model  that  can  reasonably  reproduce  the  simulation  results  for  a 
set  of  inputs,  but  at  a  fraction  of  the  computational  cost.  This  simplified  model  is  then 
coupled  to  (i.e.,  used  by)  an  optimizer  to  carry  out  the  design  optimization.  The 
optimized  design  can  later  be  verified  by  executing  the  complex  simulation  using  the 
optimal  design  variables.  The  approximation  assisted  optimization  process  is  shown 
in  Figure  1.1b.  Step-1  is  the  design  of  experiments  process,  Step-2  involves  the 
evaluation  of  the  CFD  tool  for  the  samples  generated  in  the  DOE  phase,  Step-3  is  the 
development  of  metamodel  and  Step-4  is  the  use  of  metamodel  for  optimization  and 
other  analysis. 
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Figure  1.1  (a)  Conventional  optimization,  and  (b)  Approximation  assisted  optimization 


Several  DOE  techniques  have  been  introduced  in  the  literature  (Myers  and 
Montgomery,  2002;  McKay  et  al.,  1979;  Currin  et  al.  1991),  such  as  classical  designs, 
space-fdling  designs,  etc.  Examples  of  metamodeling  techniques  include  response- 
surface  based  models  (Myers  and  Montgomery,  2002),  statistical  methods  (Cressie, 
1993),  radial  basis  functions  (Dyn  et  al.,  1986)  and  artificial  neural  networks.  Since 
the  metamodel  is  developed  with  the  (input,  output)  tuple  generated  by  the  simulation 
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in  the  DOE  phase,  clearly  this  DOE  phase  has  a  profound  effect  on  the  success  or 
failure  of  a  metamodel.  Moreover,  the  choice  of  a  DOE  method  also  depends  upon 
the  metamodeling  technique  which  will  be  used.  For  example,  classical  designs  are 
not  suitable  for  computer  simulations  and  can  sometimes  result  in  erroneous  models 
(Sacks  et  ah,  1989;  Jin  et  ah,  2001). 

Many  of  the  DOE  techniques  in  the  literature  are  such  that  the  design  points  are 
generated  in  a  single  step,  without  using  any  infonnation  from  the  corresponding 
response  or  the  resulting  metamodel.  Figure  1.2a  shows  this  non-adaptive  approach. 
Subsequently  researchers  developed  adaptive  DOE  techniques  in  which  the  DOE 
process  was  carried  out  in  stages  i.e.,  one  or  more  points  were  sampled  in  each  stage 
of  the  DOE  process.  In  this  way,  at  each  stage  a  metamodel  is  developed  and  its 
performance  is  evaluated  based  on  relevant  criteria  and  this  information  is  then  used 
to  select  points  in  the  next  stage.  It  has  been  shown  that  in  general  the  adaptive  DOE 
methods  are  superior  (see,  e.g.,  Jin  et  ah,  2002;  Li  and  Azarm,  2006;  Lam  and  Notz, 
2007)  to  non-adaptive  methods.  Figure  1.2b  shows  a  general  flow  chart  for  an 
adaptive  DOE  procedure.  As  shown  in  Figure  1.2b,  the  perfonnance  of  the 
metamodel  is  evaluated  at  each  step  (or  iteration)  and  new  points  are  sampled  based 
appropriate  performance  criteria. 
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(a)  (b) 

Figure  1.2  Generic  flow  chart  for  (a)  Non-adaptive  DOE  and  (b)  Adaptive  DOE 

An  important  aspect  of  a  computer  simulation  is  that  it  can  have  multiple  outputs. 
Thus  any  good  approximation  method  should  take  into  account  the  multiple  responses 
as  well.  One  could  conduct  a  DOE  and  develop  a  metamodel  for  each  individual 
response,  but  this  would  result  in  excessive  simulation  runs,  probably  more  than  what 
the  designer  wants  to  spare. 

A  number  of  adaptive  DOE  techniques  (Jin  et  ah,  2002,  Li  and  Azann,  2006; 
Kleijnen,  2004;  Li  et  ah,  2006)  focus  on  allocating  new  points  in  the  input  (i.e. 
design)  space  based  on  a  space  filling  criteria,  while  others  focus  on  metamodel 
performance  metric  such  as  cross-validation  only.  In  tenns  of  approximation  for 
multiple-responses,  generally,  DOE  is  seen  as  an  extension  of  a  multiresponse 
metamodel  (Romero  et  ah,  2006,  Li  et  ah,  2006)  via  the  use  of  correlation  between 
responses.  This  is  in  contrast  to  considering  the  multiple  responses  individually  (i.e., 
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developing  individual  metamodels  for  each  response)  in  the  DOE  phase.  All  of  these 
DOE  techniques,  except  for  Jin  et  al.,  (2002)  have  not  been  applied  to  design 
problems  with  higher  dimensions  (e.g.  4-8  dimensions)  nor  have  been  applied  to  a 
wide  range  of  test  problems  to  demonstrate  their  scalability.  In  some  cases  (Busby  et 
al.,  2007)  the  resulting  statistical  metamodels  have  been  compared  with  response 
surface  methods,  which  is  not  a  fair  comparison.  In  general,  the  metamodel 
performance  (i.e.,  how  well  the  metamodel  represents  the  underlying  true  function  or 
response)  can  be  improved  when  the  points  sampled  in  the  DOE  phase  systematically 
cover  the  entire  region  of  the  design  space  while  simultaneously  allocating  more 
points  in  regions  that  are  nonlinear.  Moreover,  in  multiresponse  simulations,  a  good 
DOE  technique  should  objectively  make  a  tradeoff  between  multiple  responses  as 
well  as  the  space-filling  criteria  and  the  metamodel  perfonnance  and  should  be 
scalable  to  higher  than  conventional  input  dimensions  and  to  a  wide  range  of 
problems. 

1.2  Dissertation  Objective 

The  objective  of  this  dissertation  is  (a)  to  develop  a  single-response  adaptive 
design  of  experiments  technique  for  detenninistic  black-box  computer  simulations 
that  is  scalable  and  widely  applicable,  and  an  empirical  framework  for  implementing 
the  proposed  approach  when  there  is  a  limit  on  the  computational  resources  (i.e.,  the 
number  of  simulation  calls),  (b)  to  extend  the  single-response  DOE  approach  to 
multiresponse  adaptive  design  of  experiments  technique  for  detenninistic  black-box 
computer  simulations,  and  (c)  to  develop  an  approximation-based  framework  for  the 
performance  evaluation  and  design  optimization  of  novel  air-cooled  heat  exchangers. 
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1.3  Assumptions 


The  following  assumptions  are  made  in  the  development  of  the  approach: 

(a)  The  simulation  models  are  deterministic,  i.e.,  no  matter  how  many  times 
the  simulation  is  invoked,  it  will  always  yield  the  same  output  for  the  same 
inputs. 

(b)  The  simulation  responses  are  continuous  and  the  models  are  black-box.  In 
other  words,  no  information  about  the  simulation  model  or  its  output  is 
known  before  hand  and  no  other  assumptions  are  made  with  regards  to  its 
behavior. 

(c)  The  metamodel  developed  based  on  the  proposed  approach  is  desired  to  be 
globally  (i.e.  in  the  entire  input  space)  applicable  based  on  the  available 
computational  resources. 

(d)  The  computational  resources  available  to  execute  the  simulation  are 
limited.  In  other  words,  i.e.,  the  number  of  available  simulation  calls  is 
fixed  and  the  designer  wishes  to  obtain  the  best  possible  metamodel  given 
the  DOE  points. 

(e)  The  computational  effort  for  performing  a  single  simulation  is  much 
higher  than  that  required  for  building  the  DOE  and  the  metamodel. 

1.4  Research  Thrusts 

This  dissertation  is  divided  into  three  research  thrusts  as  discussed  in  the 
following  sections. 
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1.4.1  Research  Thrust-1:  Single  Response  Adaptive  DOE 

In  this  research  thrust,  an  adaptive  DOE  approach  is  developed  for  single 

response  deterministic  computer  simulations.  An  adaptive  DOE  method  has  to  make 
tradeoffs  between  allocating  new  points  either  in  regions  of  the  design  space  that  are 
sparsely  sampled  or  in  regions  that  correspond  to  highly  sensitive  behavior  (non¬ 
linear  and  multi-modal)  of  the  response.  In  this  thrust,  it  is  shown  how  important  it  is 
to  make  a  judicious  tradeoff  between  a  space-filling  criterion  and  metamodel 
performance  in  selecting  new  experiments.  A  multiobjective  optimization  problem  is 
formulated  to  solve  this  adaptive  DOE  problem.  The  newly  developed  approach  is 
applied  to  a  suite  of  24  test  problems  to  demonstrate  its  applicability  to  a  wide  range 
of  design  problems  and  also  its  scalability  to  higher  (more  than  2)  dimensions.  The 
problems  in  the  test  suite  are  compiled  from  the  literature  and  are  challenging  (Jin  et 
ah,  2002;  Li  et  ah,  2006;  Busby  et  al.,  2007)  problems  for  adaptive  DOE  methods. 

1.4.2  Research  Thrust-2:  Multiresponse  Adaptive  DOE 

In  this  research  thrust,  an  adaptive  DOE  approach  is  developed  for  multiresponse 

deterministic  computer  simulations.  This  is  an  extension  of  the  approach  developed  in 
Research  Thrust- 1,  wherein  the  problem  is  viewed  as  a  multi-objective  optimization 
problem  of  choosing  samples  while  weighing  the  relative  perfonnance  of  the 
metamodels  for  individual  responses.  This  new  multiresponse  approach  is  applied  to  a 
total  of  10  numerical  examples  from  the  literature  and  the  results  are  compared  with 
published  results  from  previous  techniques. 
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1.4.3  Research  Thrust-3:  Approximation  Based  Framework  for  Design 
Optimization  of  Novel  Heat  Exchangers 

In  this  research  thrust,  a  new  approach  for  design  optimization  and  performance 
evaluation,  i.e.,  heat  transfer  characteristics  and  pressure  drop,  of  novel  air-cooled 
heat  exchangers  in  a  computationally  efficient  manner  is  developed.  Conventionally, 
3D  CFD  is  used  to  evaluate  the  performance  of  a  novel  air-cooled  heat  exchanger 
such  as  tube-fin  or  microchannels.  The  novel  aspect  generally  refers  to  a  new  fin 
design  or  a  tube  bundle  design.  Due  to  the  computational  burden  associated  with  3D 
CFD,  a  small  section  (i.e.,  very  short  length)  of  the  heat  exchanger  is  simulated  and 
the  results  are  assumed  to  be  valid  for  the  entire  heat  exchanger.  Even  while 
simulating  a  small  section  of  the  heat  exchanger,  the  computation  time  required  is  in 
the  range  of  few  hours  to  a  few  days.  A  new  approach  is  proposed,  in  which  2D  CFD 
(instead  of  3D  CFD)  coupled  with  a  s-NTU  heat  exchanger  model  (Jiang  et  ah,  2006) 
is  used  to  evaluate  the  performance  of  heat  exchangers  in  a  computationally  efficient 
manner.  This  new  approach  is  further  enhanced  with  the  use  of  the  techniques 
developed  in  Research  Thrusts- 1  and  2  and  the  resulting  metamodels  are  combined 
with  an  optimizer  to  find  optimal  heat  exchanger  design  configurations. 

1.5  Organization  of  the  Dissertation 

The  dissertation  is  organized  as  shown  in  Figure  1.3.  Chapter  2  provides 
definitions  and  tenninology  used  in  this  dissertation.  Also  included  in  Chapter  2  is 
background  information  on  multiobjective  optimization,  approximation  assisted 
optimization,  and  some  details  on  metamodeling  and  DOE  techniques  from  the 
literature.  Chapter  3  presents  the  Research  Thrust- 1,  wherein  a  new  approach  for 
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single  response  adaptive  design  of  experiments  is  proposed.  Chapter-4  presents  a  new 
adaptive  DOE  approach  for  multiresponse  simulations.  Chapter-5  presents  an 
approximation  based  framework  for  perfonnance  evaluation  and  optimization  of  air¬ 
cooled  heat  exchangers.  Chapter-6  concludes  the  dissertation  with  summary, 
conclusions,  research  contributions  and  directions  for  future  research. 
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Figure  1.3  Organization  of  Dissertation 
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Chapter  2  Definition  and  Terminology 


2.1  Introduction 

In  this  chapter  we  present  several  definitions  that  are  used  throughout  the 
dissertation.  We  begin  with  the  formulation  of  a  general  multiobjective  optimization 
problem  followed  by  a  discussion  of  the  resulting  set  of  solutions  and  methods  to 
obtain  solutions  to  a  multiobjective  optimization.  This  is  followed  by  an  introduction 
to  approximation  assisted  optimization  and  a  short  description  of  relevant  design  of 
experiments  and  metamodeling  approaches. 


2.2  Multi-obi ective  Optimization 


A  multi-objective  optimization  problem  (MOOP)  can  be  fonnulated  as  follows: 


minimize 

X 

ft(x) 

7  =  1,.. 

subject  to: 

gj(x)<0 

j  =c 

..,J 

4  ^  xk  ^  4 

k  =  1,. 

,.,d 

(2.1) 


where  f  refers  to  the  7th  objective  of  the  problem,  x  is  a  vector  of  d  design  variables, 
i.e.,  x  =  (xi,x2,. . .rr(i)7,  where  T  indicates  a  transpose.  There  are  M  objectives  that  are 
to  be  minimized.  It  is  assumed  that  at  least  two  of  the  M  objectives  are  conflicting  in 
nature.  The  functions  are  J  inequality  constraints.  The  vectors  x'k  and 

xuk  indicate  the  lower  and  the  upper  bounds  of  the  variable  x*. 
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2.2.1  Simulation 

A  computer  program  that  can  calculate  the  values  of  f  and  gj  is  termed  as  a 
simulation.  In  this  dissertation  we  focus  on  simulations  that  are  reproducible  or 
deterministic.  By  deterministic  we  mean  that  no  matter  how  many  times  the 
simulation  is  invoked,  for  the  same  value  of  input  design  variables,  it  always  yields 
the  same  outputs  (i.e.,  values  for  f  and  gj).  A  simulation  can  have  a  closed  (explicit) 
form  or  can  be  iterative  in  nature.  In  this  dissertation  the  simulation  is  treated  as  a 
black-box  object  as  shown  in  Figure  2.1,  where  an  example  of  single  response  (or 
output)  simulation,  Figure  2.1a,  and  multiresponse  simulation,  Figure  2.1b  is 
provided. 


Xi 


— ► 

x2 

xd 

— ► 

Black-Box 
Simulation 
(Single  Response) 


y  l 


Xi 


x2 


xd 


Black-Box 

Simulation 

(Multi-Response) 


(a)  (b) 

Figure  2.1  Black-box  simulations,  (a)  Single  response  and  (b)  Multiresponse 
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2.2.2  Input  Space,  Experiment  and  Design 

The  J-dimensional  space  formed  by  all  possible  values  of  x  e  Rd  is  termed  as  the 

design  variable  space  or  design  space.  Since  x  is  an  input  to  the  simulation,  its 
corresponding  space  is  termed  as  the  input  space.  Each  element  of  x  in  the  input  space 
is  assumed  to  be  bounded  by  a  lower  and  an  upper  bound  (see  Eq.  2.1)  and  thus  can 
be  linearly  normalized  in  the  range  [0,1],  The  output  from  the  simulation  is  termed  as 
a  true  response  when  the  full  (or  original  model)  for  the  simulation  is  used.  If  another 
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technique  such  as  a  metamodel  is  used  to  obtain  the  response  for  a  simulation,  the 
corresponding  response  is  termed  as  a  predicted  response. 

A  point  x;,(  a  specified  vector  of  x\ ,...jCd)  is  referred  to  as  an  experiment  if  it  has 
been  evaluated  (or  observed)  for  its  true  response,  i.e.,  by  invoking  the  simulation.  In 
contrast,  a  point  that  has  not  been  evaluated  for  its  true  response  is  termed  as  an 
unobserved  point. 

A  set  of  experiments  is  termed  as  a  design,  denoted  by  D. 

The  process  of  choosing  the  various  experiments  is  termed  as  sampling. 

2.3  Dominance  and  Pareto  Set 

In  multiobjective  optimization,  when  there  are  two  or  more  conflicting  objectives, 
there  cannot  be  a  single  optimum  solution  which  simultaneously  optimizes  all 
objectives.  Thus  the  solution  obtained  from  a  multiobjective  optimization  problem  is 
a  set  of  optimal  solutions  with  varying  degrees  of  objective  values. 

Considering  Eq.  (2.1),  a  solution  x\  is  said  to  dominate  (Goldberg,  1989;  Deb, 
2001)  a  solutions  xi,  when  both  (a)  and  (b)  below  hold: 

(a)  Solution  x\  is  better  than  or  equal  to  x^  in  tenns  of  all  the  objectives  and 

(b)  Solution  x\  is  strictly  better  than  X2  in  at  least  one  of  the  objectives. 

When  comparing  two  solutions,  when  the  first  condition  is  not  satisfied,  the  two 
solutions  are  said  to  be  non-dominated  with  respect  to  each  other.  In  other  words,  if  a 
point  is  not  dominated  by  any  other  point,  it  is  said  to  be  non-dominated. 
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Amongst  a  set  of  solutions  P,  the  non-dominated  set  of  solutions  P’  are  those  that 
are  not  dominated  by  any  other  point  in  P.  When  P  is  the  entire  search  space,  then  the 
resulting  non-dominated  set  P’  is  termed  as  the  Pareto  optimal  set  and  the  solutions 
are  said  to  form  a  Pareto  frontier  in  the  objective  space  (i.e.,  /-space).  Several 
algorithms  are  presented  in  Deb  (2001)  to  efficiently  identify  a  non-dominated  set 
amongst  a  given  set  of  solutions. 


§ 

- 


i4 


f 


Dominated 

Points 


Pareto 

Points 


f  • 


(minimize) 

Figure  2.2  Dominated  points  and  Pareto  points 


Note  that  in  MOOP’s  where  the  objectives  are  not  conflicting,  the  resulting  Pareto 
optimal  solution  set  will  contain  only  one  optimal  solution  (if  it  is  unique). 

2.4  Methods  for  solving  Multi-objective  optimization  problems 

2.4.1  Overview 

In  general  there  are  two  classes  of  methods  for  solving  multi-objective 
optimization  problems.  They  are:  (a)  classical  methods  and  (b)  non-classical  methods. 
Classical  methods  are  generally  gradient-based  or  direct  search  methods.  These 
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methods  are:  weighted-sum  method  (Cohon,  1978),  s-constraint  method  (Haimes, 
1971),  weighted  metric  method  (Miettinen,  1999),  value  function  method  (Keeny  and 
Raiffa,  1976;  Miettinen,  1999),  Schaffler’s  stochastic  method  (Schaffler  et  ah,  2002; 
Shukla  et  al.,  2005),  normal  boundary  intersection  method  (Das  and  Dennis,  1998), 
goal  programming  (Charnes  et  ah,  1955;  Ignizio,  1976;  Steuer,  1986)  etc.  Gradient 
based  methods  are  deterministic  in  nature  and  yield  Pareto  solutions  one  point  at  a 
time. 

The  second  class  of  methods  is  non-classical  based  methods  for  solving  MOOP’s. 
Many  of  the  non-classical  methods  are  based  on  natural  phenomenon.  These  methods 
are  generally  population  based  such  as  evolutionary  algorithms  (Goldberg,  1989; 
Deb,  2001;  Coello  Coello  et  al.,  2007),  particle  swann  optimizers  (Coello  Coello  et 
al.,  2004),  multiobjective  simulated  annealing  (Serafini,  1992;  Nam  and  Park,  2000) 
and  TimmeTs  population  based  method  (Timmel,  1980;  Shukla  et  al.,  2005)  to  name 
a  few. 

Several  variations  of  population  based  multiobjective  optimization  evolutionary 
algorithms  have  been  reported  in  the  literature  (Fonseca  and  Fleming,  1993,  1998a, 
1998b;  Srinivas  and  Deb,  1994;  Horn  et  al.,  1994;  Zitzler  and  Thiele,  1998;  Deb  et 
al.,  2000,  Deb,  2001;  Coello  Coello  et  al.,  2007).  These  methods  try  to  assign  a 
goodness  or  fitness  to  a  design  point,  based  on  its  objective  and  constraint  values.  It  is 
also  important  to  note  that  population-based  methods  require  numerous  function  calls, 
at  times  several  thousands,  of  the  simulation  to  evaluate  the  objectives  and 
constraints.  In  spite  of  the  recent  advances  in  the  field  of  high  perfonnance  and  high 
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throughput  computing,  the  use  of  population-based  optimization  methods  for  certain 
complex  problems  is  not  feasible. 

In  this  dissertation  a  multiobjective  genetic  algorithm  (MOGA)  is  used  for  solving 
MOOP’s.  It  is  important  to  point  out  that  even  though  MOGA  is  used  for  solving 
optimization  problems  in  this  dissertation,  the  approaches  proposed  are  not  limited  to 
MOGA  and  any  other  suitable  multiobjective  optimization  technique  can  be  used. 
The  details  of  a  genetic  algorithm  (GA)  and  a  MOGA  are  discussed  in  the  following 
section. 

2.4.2  Genetic  Algorithms 

Genetic  Algorithms  (GA)  as  defined  by  Goldberg  (1989)  are:  “search  algorithms 
based  on  natural  selection  and  natural  genetics”.  GA  maintains  a  pool  of  candidate 
solutions,  each  of  which  is  assigned  a  fitness  based  on  its  usefulness  or  ‘payoff. 
Fitness  is  a  scalar  measure  of  how  well  a  particular  candidate  solution  satisfies  the 
given  problem  objectives.  At  each  iteration  or  generation  of  GA,  candidate  solutions 
are  selected  for  reproduction  based  on  their  fitness  to  form  new  offspring  or  solutions. 
The  reproduction  process  is  carried  out  via  the  use  of  genetic  operators  such  as 
selection,  crossover  and  mutation  (Goldberg,  1989;  Deb,  2001).  A  set  of  probabilistic 
rules  determines  whether  a  candidate  solution  undergoes  crossover  or  mutation  and  at 
what  point.  A  powerful  feature  of  GA’s  is  that  they  search  in  multiple  directions 
simultaneously  and  do  not  require  any  gradient  information  or  other  supplementary 
information  (convex  objective  or  linear  constraints)  about  the  problem  at  hand,  only 
the  usefulness  or  the  payoff  value.  This  makes  the  GA  an  ideal  tool  for  optimization 
of  highly  non-linear  functions  involving  a  combination  of  continuous  and  discrete 
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design  variables  and  black-box  functions.  A  schematic  of  the  GA  design 
representation  and  operators  is  shown  in  Figure  2.3  for  a  heat  exchanger  optimization 
(Aute  et  ah,  2004)  example.  There  are  three  design  variables  viz.  number  of  tubes,  Nt 
(discrete),  tube  length,  L  (continuous)  and  fins  per  inch,  FPI  (discrete).  A 
chromosome  representing  one  design  is  shown  in  Figure  2.3a.  The  genetic  algorithm 
operators  of  crossover  and  mutation  are  shown  in  Figure  2.3b. 
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Figure  2.3  Genetic  Algorithms,  (a)  design  representations  and  (b)  operators 


2.4.3  Multiobjective  Genetic  Algorithm 

Since  GA’s  deal  with  a  scalar  fitness  value,  usually  some  kind  of  aggregating 

method  or  a  utility  function  is  used  to  assign  fitness  to  solutions.  A  genetic  algorithm 
that  is  modified  to  include  multiple  objectives  is  termed  as  a  Multiobjective  Genetic 
Algorithm  (MOGA).  The  main  difference  between  MOGA  and  GA  is  the 
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introduction  of  the  fitness  assignment  scheme  that  evaluates  the  goodness  of 
candidate  solutions  based  their  respective  objective  and  constraint  values.  Typically 
MOGA  considers  three  criteria  in  fitness  assignment,  viz.,  (a)  objective  values,  (b) 
constraint  values  and  (c)  crowding.  The  objective  values  are  handled  using  some  kind 
of  non-dominated  sorting  procedure  and  address  the  aspect  of  reaching  as  close  as 
possible  to  the  true  Pareto  frontier.  The  constraint  violations  are  handled  using  a 
penalty  or  custom  selection  operator  and  address  the  aspect  of  feasibility.  The 
crowding  metric  is  handled  using  some  form  of  distance  measure  such  as  niche 
(Goldberg,  1989)  or  crowding  distance  (Deb,  2000)  and  avoids  clustering  of 
solutions.  Figure  2.4  shows  a  flowchart  for  a  simple  multiobjective  genetic  algorithm. 
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Figure  2.4  MOGA  flow  chart 


2.5  Approximation  Assisted  Optimization 


2.5.1  Overview 

As  mentioned  in  Section  2.4.1,  optimization  requires  multiple  invocations  of  the 
simulation  tool  to  calculate  the  objectives  and  the  constraints.  Thus  optimization  may 
not  be  feasible  when  the  simulation  is  very  complex  such  as  a  CFD,  which  can 
execute  for  several  days  to  evaluate  one  set  of  inputs.  In  order  to  reduce  the 


20 


computational  burden  associated  with  optimization,  approximation-assisted 
optimization  is  proposed  by  researchers  (see,  e.g.,  Simpson  et  al.,  2001).  In 
approximation-assisted  optimization,  a  simplified  model  capable  of  reasonably 
representing  the  simulation  behavior  is  used  along  with  the  optimizer  to  carry  out  the 
optimization.  This  simplified  model  is  termed  as  a  metamodel  or  a  surrogate  model. 
As  explained  in  Figure  1.1  earlier,  approximation-assisted  analysis  and  optimization 
in  engineering  design  involves  four  steps:  (1)  a  design  of  experiments  (DOE)  or 
sampling  stage  where  a  design  is  chosen  for  conducting  an  experiment  to  generate 
response  data,  (2)  metamodeling  stage  in  which  the  data  from  first  step  is  used  to 
build  a  metamodel,  (3)  verification  of  the  metamodel  to  assess  how  good  the 
metamodel  is,  and  (4)  optimization  using  the  metamodel.  The  metamodel  that  is 
developed  is  computationally  more  efficient  (less  expensive)  than  the  actual 
simulation  which  it  will  replace.  An  important  step  in  approximation-assisted 
optimization  is  the  verification  of  the  optimum  solutions,  i.e.,  comparison  of  the 
metamodel  predictions  with  the  true  responses  for  the  optimum  designs.  Once  a 
verified  metamodel  is  developed,  it  can  be  used  for  multiple  optimization  studies. 

In  the  metamodel  development  phase,  the  simulation  is  invoked  multiple  times  for 
carefully  chosen  points  from  the  input  space.  This  process  is  termed  as  the  design  of 
experiments  or  DOE.  These  inputs  along  with  their  true  responses  are  then  used  to 
develop  a  metamodel.  The  design  of  experiments  is  the  structured  and  organized 
method  for  detennining  the  relationship  between  factors  or  inputs  (to  a  simulation) 
and  the  output  of  the  simulation.  The  process  of  choosing  the  various  experiments  is 
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termed  as  sampling.  Different  sampling  techniques  have  also  been  studied  and 
presented  in  the  literature  (Jin  et  al.,  2001). 

The  second  step  in  approximation  is  to  fit  the  data  to  a  model.  These  techniques 
involve  the  use  of  response  surfaces  or  other  sophisticated  stochastic  models  to 
approximate  the  true  response.  Several  computer  based  meta-modeling  approaches 
are  discussed  in  Simpson  et  al.  (2001)  and  Wang  and  Shan  (2007).  The  different 
approaches  can  be  broadly  categorized  as  (a)  response  surface  based  models  (b) 
artificial  neural  networks  (c)  statistics  based  models  such  as  Kriging  (d)  machine 
learning  and  (e)  hybrid  models.  There  are  several  metamodel  verification  methods 
proposed  in  the  literature  such  as  random-sample  evaluation,  cross-validation 
(Cressie,  1993;  Martin  and  Simpson,  2005)  etc.  The  random  sample  evaluation 
method  is  the  most  popular  and  easiest  to  implement  though  it  requires  additional 
function  evaluations.  In  this  method,  a  random  sample  of  designs  is  generated  and  the 
predicted  response  from  the  metamodel  is  compared  with  the  true  response.  The 
results  of  this  comparison  can  be  quantified  using  statistical  metrics  which  are 
discussed  in  Section  2.5.3. 

Approximation  can  be  carried  out  in  an  offline  or  online  manner.  The  difference 
between  offline  and  online  approximation  is  shown  in  Figure  2.5.  In  offline 
approximation,  as  shown  in  Figure  2.5a.,  the  design  of  experiments  and  the 
metamodel  development  is  carried  out  before  any  optimization  is  conducted.  In  online 
approximation,  shown  in  Figure  2.5b,  an  existing  metamodel  is  used  during 
optimization,  but  based  on  the  progress  (i.e.,  design  space  explored  by  the  optimizer) 
of  the  optimizer,  new  points  are  sampled  in  intermediate  stages  (e.g.,  between 
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iterations  of  the  optimizer)  and  the  metamodel  is  updated.  Online  approximation  has 
the  advantage  that  as  the  optimization  progresses  towards  optimal  solutions,  the 
accuracy  of  the  metamodels  can  be  selectively  improved  in  the  region  of  the  optima. 


(a)  (b) 

Figure  2.5  Approximation,  (a)  offline  and  (b)  online 

2.5.2  Maximum  Entropy  Design 

The  different  sampling  methods  can  be  broadly  grouped  (Simpson,  2001)  as: 
classical  methods,  space  filling  methods  and  sequential  or  adaptive  methods.  One 
such  relevant  spacing  filling  method  is  Maximum  Entropy  Design  (MED)  method 
and  is  discussed  in  this  section. 

Entropy  as  a  measure  of  information  content  was  first  introduced  by  Shannon 
(Shannon,  1948).  Lindley  (1956)  interpreted  Shannon’s  entropy  H  as  the  amount  of 
information  retrieved  from  an  experiment.  It  was  then  showed  (Shewry  and  Wynn, 
1987;  Currin  et  ah,  1988)  that  the  entropy  criterion  H,  selects  a  new  experiment  as  the 
one  that  maximizes  the  expected  retrieved  information  due  to  the  new  experiment. 
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Mathematically,  given  a  current  design  (jci ,X2,...,Xn),  the  new  experiment  is  chosen  as 
follows: 

xn+x  =argma  xH(x1,x2,...,xn;x)  (2.2) 

where  (xi,X2,...,xn)  are  the  existing  n  experiments  and  the  “argmax”  denotes  the 
optimal  solution  xn+i  of  the  maximum  entropy  optimization  problem. 

Further,  under  the  assumption  of  normal  priors  (Shewry  and  Wynn,  1987; 
Koehler  and  Owen,  1996),  it  was  shown  that  the  entropy  criterion  is  the  same  as 
maximizing  the  determinant  of  the  prior  covariance  matrix  R ,  i.e., 

xn+l  =  arg  max  det(i?)  (2.3) 

where  det  indicates  a  determinant  and  R  is  the  ((n+1)  x  (n+1))  covariance  matrix  of  x. 
Each  element  of  R  is  calculated  using  the  augmented  design  (xi,X2,....,xn,Xn+i),  where 
there  are  n  existing  designs  and  xn+i  is  the  new  candidate  design.  The  details  of 
covariance  matrix  based  on  the  normal  priors  are  given  in  Section  2.5.3. 


(a) 


(b) 


Figure  2.6  (a)  Initial  designs  and  (b)  maximum  entropy  designs 
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Figure  2.6a  shows  an  initial  design  with  5  points  in  a  two  dimensional  input  space 
and  Figure  2.6b  shows  15  additional  designs  obtained  using  the  maximum  entropy 
criterion.  Note  that  a  significant  number  of  designs  are  on  the  boundaries  of  the  input 
space  (Currin  et  ah,  1991). 

2.5.3  Kriging  Method 

Kriging  is  a  widely  used  spatial  correlation  and  prediction  technique  in  the  field 
of  Geostatistics  (Cressie,  1993;  Armstrong  1998)  and  is  named  after  the  South 
African  mining  engineer  D.  G.  Krige.  The  basic  idea  behind  Kriging  is  to  use  existing 
responses  at  some  points  x-u  where  /=1,. to  predict  the  response  at  some  unknown 
point  X(). 

Suppose  we  want  to  predict  the  model  response  at  some  point  xo.  Initially,  when 
there  are  no  samples,  then  there  is  some  assumed  uncertainty  in  this  function  value. 
This  uncertainty  can  be  modeled  as  a  random  variable  Y(x)  with  a  mean  /i  and  an 
associated  variance  a".  In  other  words,  the  function  has  a  typical  value  of  y  but  has 
some  local  variations  in  a  certain  range  as  we  traverse  the  x-space.  Consider  two 
points  Xi  and  jq.  Assuming  that  the  function  being  modeled  is  continuous  then  the 
corresponding  function  values  y(x j)  and  y(x j)  will  be  close  if  the  distance  between  x, 
andxj  i.e.  ||xi-Xj|j  is  small.  The  Kriging  approach  approximates  a  true  model  response 
v(x)  with  a  stochastic  process  as  follows: 

y(x)  =  f(x)  +  Y(x)  (2.4) 

where  y(x)  is  the  unknown  function  that  is  being  modeled  and  Y(x)  is  a  normally 
distributed  Gaussian  process.  Compared  to  ordinary  a  least  squares  method,  where  we 
are  interested  in  finding  the  functional  form  which  best  represents  the  observed  data, 
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in  Kriging,  we  are  interested  in  predicting  how  the  function  typically  behaves  and  not 
the  exact  functional  fonn.  Various  functional  forms  for  /(x)  and  Y(x)  are  available  in 
the  literature  (Jones,  2001;  Simpson  et  ah,  2001;  Martin  and  Simpson,  2005).  The 
term  f{x)  is  similar  to  a  polynomial  model  in  a  response  surface  method  and  is 
equivalent  to  a  global  mean  p  for  the  model.  The  global  mean  p  is  the  mean  of  all  the 
responses  in  the  current  design.  For  simplification  purposes  y(x)  is  generally  treated  as 
a  constant.  The  7(x)  term  represents  the  local  deviations  from  the  global  mean 
obtained  by  interpolating  the  available  data  based  on  distance  between  the 
unobserved  point  xo  and  the  sampled  points. 

The  tenn  Y(x)  is  represented  through  the  use  of  one  of  many  correlation  functions. 
One  of  the  widely  used  correlation  functions  (Sacks  et  ah,  1989;  Jones,  2001)  is: 

CorriYix^JiXj)]  =  exp(-^6>  |  xu  -xjt  \p ')  (2.5) 

/= 1 

In  this  equation,  d  is  the  dimension  of  vector  x,  xu  and  xp  are  the  /th  components  of 
the  vectors  x,  and  xh  0\  is  the  degree  of  correlation  between  the  responses  in  the  Ith 
coordinate  and  is  tenned  as  the  correlation  parameter  in  the  /th  direction,  and  p 
controls  the  smoothness  of  the  function  in  the  Ith  direction.  The  terms  0\  and  pi  provide 
a  means  for  adjusting  the  relative  importance  in  each  dimension  of  the  input  space. 
For  simplification,  a  single  value  of  6  is  used  and  the  distance  term  is  replaced  by  the 
Euclidean  distance  between  Xi  and  xy  When  one  value  of  6  is  used,  the  model  is 
tenned  as  an  isotropic  model,  which  treats  all  dimensions  equally.  In  Eq.  (2.5)  when 
p=  1,  the  correlation  is  known  as  the  exponential  correlation. 

The  effect  of  6  and p  (Jones  et  ah,  1998)  on  the  correlation  is  shown  in  Figure  2.7. 
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Figure  2.7  Correlation  vs.  Euclidean  distance  for  different  6  and  p  values 

As  observed  from  Figure  2.7,  the  correlation  tends  to  zero  as  the  distance  ||xi-Xj|| 
increases.  As  6  increases,  the  correlation  drops  rapidly  and  p  controls  the  smoothness 
of  the  correlation.  In  practice,  6  values  are  calculated  for  each  dimension  in  the  input 
space,  and  thus  when  the  correlation  is  highly  sensitive  to  one  of  the  dimensions;  the 
corresponding  6  value  will  be  higher. 

Suppose  there  are  n  responses  in  the  current  design,  then  let  y  represent  the  set  of 
n  observed  true  responses  as  given  in  Eq.  (2.6): 


Mi) 


y  = 


y(xn ) 


(2.6) 


The  uncertainty  in  the  function  values  (local  deviations)  at  the  n  points  can  be 
represented  by  a  vector  of  random  variables  Y(x)  as: 
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Y(x)  = 


(2.7) 


Y(x  i) 

Y(xn)_ 

This  vector  has  a  covariance  matrix  Cov  given  by: 


Cov(Y )  =  aR 


(2.8) 


where  R  is  an  n  x  n  correlation  matrix  with  the  (i,j)  element  given  by  Eq.  (2.5).  The 
diagonal  elements  of  R  are  always  of  the  form  Corr[Y  (x,),  Y  (x,)]  and  thus  are  always 
equal  to  1 .  Let  /  denote  an  n  x  1  vector  of  ones  and  r  denote  the  correlation  of  T(xo), 
the  unobserved  point,  with  T(x,),  the  current  designs,  as 


Corr[T(x0),T(x1)] 


Corr[Y  (x0  ),  Y  (xn  )] 


(2.9) 


The  values  of  the  correlation  parameters,  such  as  pi,  a ,  0  and  pi  need  to  be 
estimated.  They  are  obtained  by  maximizing  the  likelihood  function  or  in  other 
words,  to  model  the  functions  behavior  so  that  it  closely  represents  the  observed  data. 
Maximizing  the  likelihood  function  provides  an  estimate  of  the  optimal  values  of  pi, 
cf,  as  functions  of  R  as  follows 


fi  = 


ITR-ly 

I'R  7 


(2.10) 
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and 


.2  _  iy-ltf  'R~\y-Iju) 
n 


(2.11) 


Then  the  estimated  response  for  an  unobserved  point  x0  is  given  using  the 
Kriging  predictor  as: 


y(xQ)  =  ju  +  rTR  \y- 1 ju) 


(2.12) 


In  addition  to  the  above  predictor,  the  Kriging  measure  of  uncertainty  in  the 
estimated  response  can  be  calculated  as  follows: 


s2(x0)  =  &2 


1  -rTRlr  + 


(1  -rTR~lrf 
IR  7 


(2.13) 


The  derivation  for  this  standard  error  is  provided  in  Sacks  et  al.  (1989)  and  Jones 
(2001).  It  can  be  seen  that  the  variance  s  is  zero  for  an  observed  point.  In  utilizing 
Kriging  predictor  in  approximation  assisted  optimization  (Bakker,  2000;  Jones, 
2001),  this  standard  error  can  serve  as  a  basis  for  making  the  decision  between  using 
the  predicted  response  and  invoking  the  analyzer  functions  to  obtain  a  true  response. 
The  prediction  of  the  standard  error  is  a  big  advantage  of  Kriging  over  other 
metamodeling  methods  since  the  metamodel  can  then  be  dynamically  updated  based 
on  the  responses  during  a  given  optimization  procedure.  Furthermore,  as  mentioned, 
Kriging  does  not  require  a  functional  form,  though  the  choice  of  the  correlation 
function  is  problem  dependent.  Simpson  et  al.  (2001)  found  that  Kriging  is  extremely 
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flexible  and  suitable  for  detenninistic  computer  experiments  and  recommend  the  use 
of  Kriging  metamodels  when  the  number  of  input  variables  is  less  than  50. 


2.5.4  Metamodel  Performance  Metrics:  Errors  Verification 

The  performance  of  a  metamodel  is  generally  evaluated  by  comparing  the 

response  obtained  using  the  metamodel  with  the  true  response  for  a  given  set  of 
points  termed  as  a  test  sample.  This  test  sample  can  be  chosen  randomly  or  using 
some  systematic  method  e.g.,  based  on  a  statistical  distribution. 

Let  y(xj)  be  the  true  response  from  a  simulation  and  y(xi )  be  the  predicted 

response  for  a  test  sample  of  size  n,  where  i=l,...,n.  Then  we  can  define  a  global 
measure  of  the  accuracy  of  the  metamodel  tenned  as  the  root  mean  square  error 
(RMSE)  as  follows: 

Rmse  =  J-ily(x,)-y(x,)]2  (2.14) 

M  n  /=1 

The  maximum  absolute  error  (MAE)  represents  the  maximum  deviation  of  the 
predicted  response  from  the  true  response  for  all  the  test  samples  and  is  given  by: 

MAE  =  max(|  y(x,)~ y(xt)  \ ),i  =  1 (2.15) 

In  addition,  we  can  define  RMSE  and  MAE  based  on  relative  errors.  The  relative 
error  for  a  sample  x\  can  be  defined  as: 


_  fa) 
y(y) 


(2.16) 


Then,  the  corresponding  relative  RMSE  and  relative  MAE  can  be  defined  as: 
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RRMSE  = 


(2.17) 


and 

RMAE  =  max(|  e;.  |)  (2.1 8) 

The  RMSE  and  MAE  consider  only  the  numerical  magnitude  of  the  errors.  The 
relative  RMSE  and  MAE  metrics  are  useful  when  the  numerical  range  of  the  response 
i.e.,  difference  between  the  minimum  and  the  maximum  values  differs  by  several 
orders  of  magnitude.  RRMSE  and  RMAE  are  useful  in  practical  engineering 
examples  as  demonstrated  in  Chapter  5. 

Even  though  the  above  metrics  or  RMSE  and  MAE  are  widely  used  in  the 
literature  for  assessing  metamodels,  these  measures  are  sensitive  to  the  sample  size. 
Hamad  (2006)  introduced  a  new  metric  termed  as  metamodel  acceptability  score 
(MAS)  that  is  less  sensitive  to  the  sample  size  and  is  more  intuitive  in  the  context  of 
the  use  of  metamodels  for  global  optimization.  The  designer  needs  to  provide  a 
bound  on  the  acceptable  predictor  error  of  the  metamodel.  Let  £MAS  be  an 
acceptability  threshold,  then  a  predicted  response  y  is  said  to  be  acceptable  if  Eq. 
(2.19)  is  satisfied: 

(1-^s)  <  (t,/v;)<(1+^mas)  (2.19) 

If  there  are  n  test  samples,  of  which  m  responses  are  acceptable  as  per  above 
equation,  then  MAS  is  defined  as: 

m 

MAS  =  —  x  100  (2.20) 

n 
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The  acceptability  criteria  is  a  fraction,  e.g.,  £\,(AS  =0. 1  implies  a  tolerance  of  10% 

in  the  predicted  response.  It  should  be  pointed  out  that  the  acceptability  criteria  need 
not  be  symmetric.  This  is  important  especially  when  certain  information  about  the 
application  is  known  before  hand.  For  example,  in  an  air-cooled  heat  exchanger 
design,  a  metamodel  that  over  predicts  the  air-side  pressure  drop  (in  a  reasonable 
range)  could  still  be  acceptable,  since  this  will  result  in  a  conservative  response  and 
may  result  in  computational  savings.  A  numerical  issue  with  Eq.  (2.19)  is  the  case 
wheny;  =0.  Hamad  (2006)  proposes  the  following  remedies:  (a)  if  the  number  of 

cases  with  zero  true  response  is  a  small  fraction  of  the  overall  sample  size,  then  these 
cases  can  be  neglected,  (b)  if  the  former  is  not  the  case,  then  a  ‘shift-transformation’ 
can  be  applied  to  yn  i.e.,  y.  =yi  +  8  to  avoid  with  division  with  zero.  The  shift 
transformation  should  not  affect  the  overall  value  of  the  MAS. 

It  should  be  pointed  out  that  in  the  case  of  multiple  responses,  the  metrics  of 
MAS,  RMSE,  MAE,  RRMSE  and  RMAE  are  calculated  individually  for  each 
response. 

Another  method  of  metamodel  verification  is  the  method  of  cross-validation 
(CV)  error  calculation  (Cressie,  1993;  Meckesheimer  et  ah,  2002).  A  p-fold  cross- 
validation  approach  is  carried  out  as  follows:  (i)  the  existing  design  (i.e.,  inputs  and 
responses)  are  split  into  p  different  mutually-exclusive  subsets,  (ii)  a  metamodel  is  fit 
p  times,  each  time  leaving  out  one  subset,  the  resulting  metamodel  is  used  to  predict 
the  response  for  the  omitted  subset  and  an  error  measure  is  calculated.  This  error 
measure  is  the  cross-validation  error  (CVE)  and  is  defined  in  the  literature  as  the 
absolute  difference  between  the  predicted  and  the  actual  response.  A  special  case  of 


32 


cross-validation  method  is  the  leave-one-out  (LOO)  error  calculation  and  this  is 
discussed  in  more  detail  in  Chapter  3.  For  Kriging  metamodels  Meckesheimer  et  al. 
(2002)  propose  to  use  0.  In  or  n  as  the  value  for  p,  where  n  is  the  number  of  samples 
in  the  current  design.  We  note  that  while  0.1  n  is  a  good  number  for  assessing  the 
metamodel  performance,  it  is  computationally  very  expensive  as  n  increases,  since 
each  cycle  of  CV  calculation  requires  the  metamodel  to  be  refitted.  The  advantage  of 
this  approach  is  that  it  does  not  require  any  additional  experiments.  In  statistics  a 
technique  called  boot-strapping  (Cressie,  1993)  is  also  used  for  model  verification, 
but  in  contrast  to  cross-validation,  boot-strapping  requires  additional  function  calls. 

Throughout  this  dissertation,  the  metrics  of  RMSE,  MAE,  RRMSE,  RMAE  and 
MAS  are  used  to  evaluate  the  performance  of  metamodels. 

2.5.5  Maximin-Scaled  Distance  (MSD)  Method 

The  new  adaptive  approach  proposed  in  Chapter  3  is  compared  against  two  other 

adaptive  DOE  methods  from  the  literature.  The  first  of  these  is  the  maximin  scaled 
distance  method  proposed  by  Jin  et  al.  (2002).  The  maximin  scaled  distance  (MSD) 
approach  (Jin  et  al.,  2002)  is  an  enhanced  version  of  the  original  maximin  distance 
approach  (Johnson  et  al.,  1990)  to  make  use  of  the  information  available  from  the 
metamodel.  The  method  chooses  the  next  sample  as  the  point  that  is  farthest  from  all 
points  in  the  current  design.  The  distance  used  for  comparison  is  a  scaled  distance  as 
given  in  Eq.  (2.21).  Mathematically,  the  next  sample  xn+\  is  chosen  as  the  solution  of 
the  problem: 
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Xj^X: 


xn+l  =  arg  max  [min  ds(xt ,  x . )] 
where 


k  =  1, 


(2.21) 


The  weight  parameters  ah  are  chosen  to  be  the  same  as  the  correlation  parameters 
9,  obtained  from  a  Kriging  metamodel.  These  weight  parameters  reflect  the  relative 
importance  of  the  individual  design  variables.  At  each  stage,  a  Kriging  metamodel  is 
built  and  the  estimated  0d  values  are  used  to  solve  the  optimization  problem  from  Eq. 
(2.21)  to  obtain  the  next  sample. 


2.5.6  Accumulative  Error  (ACE)  Method 

The  second  adaptive  DOE  method  from  literature  is  the  accumulative  error  (ACE) 

method  proposed  by  Li  and  Azarm  (2006)  and  in  Li  (2007).  The  new  approach 
proposed  in  Chapter  3  was  inspired  in  part  by  the  ACE  method. 

The  ACE  method  develops  a  metamodel  for  leave  one  out  cross  validation  errors 
using  a  simple  isotropic  correlation  along  with  the  use  of  a  distance  threshold  to  avoid 
clustering  of  points.  Isotropic  implies  that  the  functional  behavior  is  the  same  along 
each  dimension  of  the  input  space.  In  the  general  the  cross-validation  error 
metamodel  tends  to  predict  higher  error  values  for  points  that  are  close  to  the  points 
in  the  current  design.  Let  D  be  the  current  design  with  n  points,  then  the  new  sample 
is  given  by: 
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(2.22) 


xH+i  =  arg  max  ^  e'wo [exp(-a  ||  x,  -x0 1|)] 

1=1 

s.t. 

ltx/-xoll^dCfit,xo0D,xieD 

where  dcnt  is  a  cluster  threshold  distance  calculated  based  on  the  points  in  the  current 
design  and  a  is  an  adaptive  correlation  parameter.  The  details  of  calculation  of  dcnl 
and  a  can  be  found  in  Li  and  Azarm  (2006)  and  Li  (2007).  It  should  be  pointed  out 
that  a  is  a  scalar  value  and  as  such  does  not  take  into  account  the  effect  of  individual 
inputs  on  the  response.  The  cluster  threshold  dcnt  is  used  to  ensure  that  new  sample 
points  are  sufficiently  away  of  existing  points. 

In  the  next  chapter,  a  cross-validation  based  adaptive  design  of  experiments 
method  for  single  response  deterministic  simulations  is  presented. 
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Chapter  3  Single  Response  Adaptive  Design  of  Experiments 


3.1  Introduction 

This  chapter  presents  a  new  adaptive  design  of  experiments  method  for  single 
response  deterministic  simulations.  As  introduced  in  Chapters  1  and  2,  any 
optimization  algorithm  requires  numerous  invocations  of  the  simulation  tool  which 
calculates  the  objective  and  constraint  functions  of  the  optimization  problem.  Often  in 
real  engineering  design  problems  these  functions  (or  their  corresponding  simulations) 
do  not  have  a  closed  form  expression,  i.e.,  they  are  iterative  in  nature  and  can  involve 
a  solution  of  linear,  nonlinear  and/or  differential  equations.  An  example  would  be  a 
complex  finite  element  or  a  computational  fluid  dynamics  simulation.  Even  when 
executed  on  the  most  sophisticated  computing  hardware,  these  simulations  may 
require  several  minutes  to  several  hours  to  execute  and  evaluate  for  single  set  of 
inputs.  A  technique  often  used  to  reduce  this  computational  burden  is  the  use  of 
offline  metamodeling.  The  goal  of  offline  metamodeling  is  to  develop  a  reasonably 
accurate  metamodel  of  the  entire  solution  space  (objectives  and  constraints)  or  a 
particular  region  of  interest  of  the  solution  space  using  as  few  experiments  as 
possible. 

Developing  a  reasonable  metamodel  with  the  fewest  experiments  possible  clearly 
emphasizes  the  importance  of  design  of  experiments  (DOE)  in  selecting  the 
experiments.  The  accuracy  of  the  metamodel  will  heavily  rely  on  the  DOE  method 
used.  In  this  chapter  an  adaptive  DOE  approach  for  multi-input  single  output  (or 
single  response)  simulations  is  presented.  We  refer  to  this  method  as  the  Space- 
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Filling  Cross  Validation  Tradeoff  (SFCVT)  method.  In  SFCVT,  the  problem  of 
adaptively  generating  new  experiments  at  a  given  point  is  viewed  as  a  tradeoff 
between  a  “space-filling”  criterion  and  a  Leave-One-Out  (LOO)  error  criterion. 

The  remainder  of  this  chapter  is  organized  as  follows:  Section  3.2  provides  an 
overview  of  the  previous  work  in  the  area  of  adaptive  DOE  methods.  Section  3.3 
explains  the  terminology  and  definitions  used  to  describe  the  proposed  approach. 
Section  3.4  describes  the  proposed  approach  in  detail.  In  Section  3.5,  the  proposed 
approach  is  applied  to  suite  of  24  numerical  test  problems  and  the  results  are 
discussed  in  detail.  The  proposed  approach  is  also  compared  with  three  other  DOE 
methods  from  the  literature.  An  empirical  procedure  to  allocate  available  function 
calls  to  the  proposed  adaptive  design  of  experiments  and  to  a  random  verification  of 
the  obtained  metamodel  is  proposed  in  Section  3.6.  The  approach  and  results  are 
summarized  in  Section  3.7. 

3.2  Related  Work 

Several  recent  publications  (Giunta  et  ah,  2003,  Wang  and  Shan,  2007;  Shao  and 
Krishnamurty,  2008)  provide  an  overview  of  DOE  methods  suitable  for  modern 
computational  engineering  design.  DOE  methods  can  be  classified  based  on  two 
different  criteria  viz.,  number  of  stages  or  iterations  and  adaptive  vs.  non-adaptive 
criteria. 

Based  on  the  number  of  stages,  DOE  methods  can  be  classified  as  single-stage  or 
multi-stage  methods.  In  single-stage  methods,  all  the  required  experiments  are 
sampled  in  a  single  run.  In  multi-stage  methods,  also  known  as  sequential  designs 
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methods;  one  or  more  points  are  sampled  at  each  stage.  Generally,  multi-stage 
methods  with  one  point  per  stage  are  shown  to  work  well  overall  (Romero  et  ah, 
2006). 

DOE  methods  can  also  be  classified  as  being  non-adaptive  or  adaptive.  In  non- 
adaptive  DOE  methods,  only  the  information  from  the  experiments  (i.e.,  design  or 
input  space)  is  used.  Several  non-adaptive  DOE  methods  can  be  found  in  the 
literature.  Classical  designs  are  used  when  physical  experiments  are  conducted.  Some 
examples  of  classical  designs  (Myers  and  Montgomery,  2002)  include  full  and 
fractional  factorial  designs,  central  composite  designs  and  Box-Behnken  designs.  The 
classical  methods  tend  to  allocate  points  on  or  near  the  boundaries  of  the  design  space 
and  leave  a  few  points  in  the  centre.  Sacks  et  al.  (1989)  and  Jin  et  al.  (2001)  conclude 
that  in  the  case  of  detenninistic  computer  simulations,  a  good  experimental  design 
tends  to  fill  the  entire  designs  space  instead  of  focusing  only  on  the  boundaries  or  at 
the  center.  Examples  of  space  filling  designs  include  Latin  Hypercubes  (McKay, 
1979)  abbreviated  as  LHC,  maximum  entropy  design  (Shewry  and  Wynn,  1987; 
Koehler  and  Owen,  1996),  mean  squared  error  (Jin  et  al.,  2002),  integrated  mean 
squared  error  (Sacks  et  al.,  1989),  maximin  distance  approach  (Johnson  et  al.,1990), 
orthogonal  arrays  (Taguchi,  1987;  Owen,  1992),  Hammersley  sequences 
(Kalagnanam  and  Diwekar,  1997)  which  are  also  known  as  low  discrepancy 
sequences  (Wang  and  Shan,  2007)  and  uniform  designs  (Fang  et  al.,  2000). 

A  DOE  approach  is  said  to  be  adaptive  when  the  information  from  the 
experiments  (inputs  and  responses)  as  well  as  information  from  the  metamodel  is 
used  in  selecting  the  next  sample.  Adaptive  approaches  are  typically  superior  to  non- 
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adaptive  approaches.  An  adaptive  approach  generally  begins  with  an  initial  design 
chosen  randomly  or  using  some  space-filling  method.  A  metamodel  is  constructed 
using  the  initial  experiments  and  then  new  samples  are  chosen  by  systematically 
evaluating  the  current  design  and  the  metamodel. 

The  adaptive  DOE  approaches  proposed  in  the  literature  can  be  classified  into  two 
categories:  (a)  DOE  for  globally  accurate  metamodels  and  (b)  DOE  for  optimization 
and/or  other  analysis.  The  DOE  for  globally  accurate  metamodels  are  such  that  they 
can  be  used  to  provide  metamodels  which  give  a  good  estimate  of  the  response  in  the 
entire  design  space.  On  the  other  hand,  the  DOE  for  the  optimization  category, 
involves  the  use  of  metamodels  but  they  also  add  new  samples  based  on  progress  of 
optimization,  i.e.,  estimation  of  an  optimum.  The  various  adaptive  DOE  approaches 
proposed  in  the  literature  use  one  or  both  of  the  two  criteria  as  a  basis  for  selection  of 
new  sample  points,  i.e.,:  (a)  estimate  of  an  optima  and  (b)  estimate  of  errors.  Methods 
which  use  an  estimate  of  the  optima  include  the  methods  by  Cox  and  John  (1997) 
termed  as  sequential  design  for  optimization  (SDO),  Watson  and  Barnes  (1995), 
Jones  et  al.  (1998),  Sasena,  2002,  Sasena  et  ah,  2000.  These  methods  try  to  locate  the 
optimum  using  the  metamodel  as  each  stage  and  sample  points  close  to  the  estimated 
optimum.  Methods  that  are  based  on  the  estimate  of  errors  include  those  by  Lin  et  al. 
(2004a,  2004b).  Lin  et  al.  (2004a,  2004b)  propose  to  use  a  supplementary  metamodel 
to  model  the  prediction  errors  of  the  metamodel  based  on  a  test  sample.  The  test 
sample  is  generated  based  on  a  maximum  entropy  principle,  i.e.,  maximizing  the 
determinant  of  a  modified  covariance  matrix.  They  used  Kriging  for  the  response 
metamodel  and  multivariate  adaptive  regression  splines  for  the  supplementary 
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metamodel  of  prediction  errors.  Their  approach  was  applied  to  only  one  test  problem 
with  one  design  variable  and  one  response.  Further  tests  are  necessary  to  evaluate  the 
general  applicability  of  their  method.  The  method  by  Jin  et  al.  (2002)  is  the  MSD 
method  and  the  one  by  Li  and  Azarm(2006)  is  the  ACE  method  described  in  Chapter 
2.  Goel  et  al.  (2006)  also  use  the  prediction  uncertainty  in  responses  along  with 
multiple  metamodels.  In  their  approach,  multiple  metamodels  are  used  to  model  the 
response  and  then  the  uncertainty  in  the  predicted  response  from  these  metamodels  is 
used  to  choose  new  samples.  The  final  predicted  response  is  a  weighted  combination 
of  the  responses  from  the  individual  metamodels  used  in  the  DOE  phase.  Their 
approach  was  presented  for  single  response  simulations. 

Cross-validation  error  based  methods  use  either  the  actual  cross-validation  error 
or  the  cross-validation  variance.  In  general,  a  metamodel  is  built  using  the  cross- 
validation  (CV)  error  or  a  cross-validation  variance  (CVV)  and  new  samples  are 
chosen  based  on  the  maximum  value  predicted  by  the  CV  or  CVV  metamodel. 
Approaches  based  on  CV  error  include  those  by  Jin  et  al.  (2002),  Busby  et  al.  (2007) 
and  the  Accumulative  Error  (ACE)  approach  Li  and  Azann  (2006)  and  Li  (2007), 
which  use  LOO  errors. 

The  method  by  Busby  et  al.  (2007)  uses  an  adaptive  gridding  algorithm  which 
divides  the  input  space  into  a  non-uniform  grid  cells.  For  each  of  these  grid  cells,  a 
CVE  is  calculated  based  on  the  points  that  already  exist  in  the  design.  In  the  cells 
with  maximum  CVE,  a  locally  optimal  point  is  sought  using  the  MED  method.  The 
grid  cells  that  do  not  have  any  points  are  assigned  an  arbitrarily  large  value  of  CVE, 
which  can  result  in  false  predictions  in  subsequent  stages.  The  perfonnance  of  this 
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method  was  compared  against  response  surface  based  metamodeling  techniques  and 
was  shown  to  be  superior.  Further  numerical  experiments  and  comparisons  using 
more  sophisticated  metamodels  are  required  to  evaluate  the  method’s  applicability. 

The  methods  by  Jin  et  al.(2002)  and  Li  and  Azann  (2006)  use  a  relatively  simple 
approach;  first  to  predict  the  cross-validation  errors  at  different  points  in  the  design 
space  and  then  new  points  are  chosen  based  on  maximum  value  of  the  predicted  CV 
error.  In  general  CV  error  based  approaches  tend  to  allocate  points  close  to  each  other 
resulting  in  clustering.  The  methods  by  Jin  et  al.  (2002)  and  Li  and  Azann  (2006)  can 
be  summarized  in  two  steps  as  follows: 

(a)  Develop  a  simplified  metamodel  for  the  LOO  prediction  errors  and  scale  the 
resulting  LOO  errors  with  the  distance  of  the  unobserved  points  from  existing 
samples. 

(b)  Choose  new  sample  as  the  one  with  the  maximum  value  of  the  LOO  error 
calculated  in  Step-a. 

Jin  et  al.  (2002)  use  an  averaging  of  LOO  prediction  errors  while  Li  and  Azann 
(2006)  use  a  weighted  combination.  The  approach  by  Li  and  Azann  (2006)  is 
described  in  more  detail  in  Chapter  2  (see  Section  2.5.6).  Methods  that  use  CVV 
include  those  by  Kleijnen  and  van  Beers  (2004)  and  Romero  (2006).  Kleijnen  and 
van  Beers  (2004)  uses  a  supplementary  metamodel  of  CVV  and  new  samples  are 
selected  based  on  the  maximum  value  of  the  predicted  CVV.  They  concluded  that 
their  approach  outperforms  the  LHC  design  and  the  MSE  method  when  applied  to  a 
single  response  test  problem.  Romero  et  al.  (2006)  extend  the  approach  proposed  by 
Kleijnen  (2004)  method  to  multiple  responses. 
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Another  method  from  the  literature  is  the  mean  squared  error  (MSE)  method 
(Jones,  2001;  Jin  et  ah,  2002)  which  chooses  the  next  sample  based  on  the  maximum 
value  of  the  estimated  mean  squared  error  in  the  response  predicted  by  an  existing 
Kriging  metamodel.  Mathematically, 

xM  =  arg  max  s2  (x)  (3.1) 

2 

where  5"  is  given  by  Equation  2.13. 

Hendrickx  et  al.  (2006)  propose  a  distributed  computing  approach  for  adaptive 
DOE  and  metamodeling,  in  which  multiple  metamodels  are  evaluated  at  each  stage 
simultaneously  and  the  best  model  is  chosen  for  the  next  step. 

Overall,  the  shortcomings  of  the  above  mentioned  approaches  can  be  summarized 
as  follows:  (a)  the  methods  use  estimated  optima  and/or  the  estimation  errors  in  the 
response,  which  can  be  misleading  especially  when  the  number  of  initial  samples  is 
small,  (b)  clustering  of  samples,  (c)  isotropic  models  fail  to  capture  the  response 
space  and  (d)  scalability  of  the  approaches  to  higher  dimensions  and  a  wide  range  of 
applications  is  not  demonstrated  (e)  comprehensive  comparison  with  other  relevant 
approaches  from  literature.  In  order  to  address  these  shortcomings  an  adaptive  DOE 
approach  for  single  response  detenninistic  simulations  is  proposed  in  this  chapter.  A 
test  suite  of  24  numerical  and  engineering  problems  from  the  literature  is  also 
compiled  so  that  the  applicability  and  scalability  of  the  proposed  approach  and  the 
previous  DOE  approaches  (i.e.,  MED  by  Shewry  and  Wynn,  1987,  MSD  by  Jin  et  al., 
2002,  and  ACE  by  Li  and  Azarm,  2006)  can  be  evaluated  over  a  range  of  input 
dimensions  and  response  behaviors.  The  objective  of  the  proposed  approach  is  to 
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obtain  a  good  estimate  of  the  entire  response  space.  There  are  many  applications 


where  a  metamodel  with  reasonable  accuracy  in  the  entire  design  space  is  desired, 
such  as  the  two  cases  in  the  following: 

(a)  Multiobjective  optimization:  In  multiobjective  optimization  problems,  the 
designer  is  not  trying  to  maximize  or  minimize  one  particular  objective 
function,  but  instead  trying  to  find  out  a  tradeoff  between  the  multiple 
objectives.  Such  a  tradeoff  generally  does  not  occur  at  the  extreme  values 
of  the  objective  functions. 

(b)  Multi-level  design:  Consider  a  multi-level  design  problem  as  shown  in 
Figure  3.1.  The  system  level  objective  /  is  a  function  of  sub-system 
objectives  f\  and  /2  and  also  additional  design  variables.  Approximation  is 
used  at  the  subsystem  levels  to  predict  the  responses  /  and  /.  In  these 
types  of  design  problems,  the  designer  is  not  interested  in  the  extreme 
points  of  the  function  being  approximated  (i.e.,  /  and  />);  instead  the 
designer  is  interested  in  finding  an  estimate  to  a  different  objective 
function  (i.e.,  /  in  this  case).  The  functions  being  approximated  may  not 
contribute  to  the  designer’s  objective  function  in  a  mono  tonic  manner 
(i.e.,  F  may  be  highly  non-linear)  and  thus  the  extreme  point  of  the 
designer’s  objective  function/ may  not  correspond  to  the  extreme  points 
of  the  functions /  and/  which  are  being  approximated.  A  design  problem 
of  such  nature  is  studied  in  Chapter  5. 
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Figure  3.1  Multi-level  design  problem 


3.3  Terminology 

This  section  introduces  the  specific  terminology  used  in  the  proposed  approach. 

3.3.1  Response  Characteristics 

The  response  or  the  output  of  an  engineering  simulation  can  have  a  Continuous 
And  Multi-Modal  (CAMM)  characteristic.  Figure  3.2  shows  examples  of  a  CAMM 
and  a  non-CAMM  region  in  the  response  space  for  a  single-input  single-output 
simulation.  The  response  function  used  to  generate  Figure  3.2  is  adapted  from  Li 
(2007). 
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Figure  3.2  Regions  of  sensitive  response 


An  example  of  response  which  is  highly  multi-modal  is  shown  in  Figure  3.3.  This 
plot  corresponds  to  the  Schwefel’s  function  (Hedar,  2005).  Any  adaptive  DOE 
approach  should  be  able  to  identify  such  CAMM  regions  in  the  response  space. 
Intuitively,  we  would  expect  the  approach  to  sample  more  points  in  the  CAMM 
region. 
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Figure  3.3  Highly  multimodal  response  function 

3.3.2  Leave-One-Out  Cross-validation  Error 

The  method  of  cross-validation  (Cressie,  1993;  Meckesheimer  et  ah,  2002)  has 

been  studied  extensively.  The  advantage  of  this  method  is  that  no  additional  function 
evaluations,  i.e.,  simulations  calls  are  required.  The  method  of  cross-validation  is 
used  to  assess  metamodel  accuracy  given  a  current  or  existing  design.  At  each  step, 
one  point  is  omitted  from  the  current  design  and  the  remaining  points  are  used  to 
build  a  metamodel.  The  resulting  metamodel  is  then  used  to  evaluate  the  response  for 
the  omitted  point  and  an  error  measure  (e.g.,  actual  error  or  absolute  error)  is 
calculated.  This  cross  validation  error  is  termed  as  the  Leave-One-Out  (LOO)  error. 
Let  there  be  n  experiments  in  a  design.  Then  the  LOO  error  for  a  point  Xi  in  the  design 
can  be  represented  as: 
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(3.2) 


eLoo(xi)=\  y,-y, 

It  is  important  to  point  out  that  the  LOO  error  is  a  not  a  sufficient  measure  of  the 
metamodel  accuracy.  It  is  actually  a  measure  of  the  sensitivity  or  insensitivity  of  the 
metamodel  to  the  lost  information  (Cressie,  1993).  Thus  an  insensitive  metamodel  is 
not  necessarily  accurate.  But  the  LOO  error  calculation  does  allow  us  to  probe  the 
response  space  for  CAMM  regions. 

An  example  of  the  LOO  error  calculation  and  the  resulting  metamodel  response  is 
shown  in  Figure  3.4  (adapted  from  Li,  2007).  The  true  response  for  the  function  is 
also  shown  in  the  figure.  The  metamodel  was  developed  with  10  sample  points.  Then 
the  LOO  error  (qloo)  was  calculated  by  omitting  the  point  x=0.25.  The  resulting 
metamodel  and  the  LOO  error  is  shown  in  the  figure.  It  is  important  to  note  how 
sensitive  the  metamodel  is  to  the  Tost  information’  due  to  omission  of  the  point 
x=0.25. 

Eq.  (3.2)  is  the  generally  used  fonn  of  the  LOO  error  calculation  in  the  literature. 
While  it  does  allow  us  to  probe  the  response  space  for  CAMM  regions  of  the 
response,  the  LOO  error  quantity  is  sensitive  to  the  magnitude  of  the  response 
surface.  This  issue  can  be  addressed  by  using  a  proper  nonnalization  of  the  response 
values  and  by  using  these  normalized  values  in  the  DOE  and  metamodel  development 
phases.  But  the  infonnation  required  for  normalization,  i.e.,  lower  and  upper  bounds 
on  the  response  is  not  generally  known  apriori.  To  alleviate  this  problem,  it  is 
proposed  to  use  the  following  definition  of  LOO  error  (i.e.,  relative  LOO  error): 

=  —  (3-3) 

yt 
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Figure  3.4  LOO  error  calculation 

This  relative  LOO  error  can  be  intuitively  much  more  practical,  since  in 
engineering  design  problems,  a  globally  accurate  metamodel  is  expected  to  predict 
the  response  as  close  as  possible,  generally  with  a  certain  error  bound  (e.g.,  with  10% 
of  true  value).  The  difference  between  the  actual  LOO  and  the  relative  LOO  is  shown 
in  Figure  3.5  for  a  test  problem  P8  from  Appendix-A.  The  design  comprised  of  20 
points  and  the  LOO  errors  were  calculated  for  these  points  using  Equations  (3.2)  and 
(3.3).  Note  that  relative  LOO  error  is  highest  for  a  point  close  to  (0.5,  0.5)  where  as 
the  actual  LOO  error  value  is  higher  for  several  points,  away  from  point  (0.5,  0.5), 
especially  on  the  boundary.  Similar  trends  are  observed  for  other  test  problems  from 
Appendix-A. 
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Thus,  the  conventional  (absolute)  LOO  error  can  be  misleading  and  its  use  as  a 


criterion  for  selecting  new  samples  can  lead  to  inefficient  designs.  In  this  dissertation 
the  proposed  relative  LOO  error  is  used  and  henceforth  the  LOO  error  refers  to  the 
relative  LOO  error  given  by  Eq.  (3.3). 


(c) 


(d) 


Figure  3.5  LOO  Error  (a)  Response  surface  for  P8,  (b)  initial  design,  (c)  Actual  LOO 

error  and  (d)  relative  LOO  error 
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3.3.3  Level  of  Dominance 

As  mentioned  in  Chapter  2,  all  the  solutions  in  a  Pareto  set  are  non-dominated 
with  respect  to  each  other.  This  notion  of  non-domination  does  not  consider  the 
magnitude  of  different  objective  values.  In  the  context  of  multi-criteria  decision 
making,  Koppen  et  al.  (2005)  introduce  the  concept  of  fuzzy  dominance  using  the 
metric  termed  as  degree  of  dominance.  The  degree  of  dominance  (DoD)  first 
introduced  in  Koppen  et  al.  (2005),  is  defined  as  follows:  Let  /(1)  and  /(2)  be  two 
objective  values  of  two  solutions  in  the  context  of  a  multi-objective  optimization 
problem  with  M  objectives.  Then  the  degree  of  dominance  of  solution  /(1)  over 

solution  /(2)  is  a  constrained  norm  in  the  objective  space  and  is  defined  as: 

M 

Drain 

«2=J=!— i -  (3.4) 

nr 

i=i 

jul2  is  the  amount  by  which  solution  (1)  dominates  solution  (2).  Eq.(3.4)  is  a 
constrained  division,  i.e.,  if  any  of  the  f-  :  is  zero,  then  the  corresponding  indices  are 
ignored  from  the  calculation.  //2I  can  be  defined  in  a  similar  manner.  The  above 
definition  has  the  following  interesting  properties: 

(a)  The  degree  of  dominance  is  not  symmetric  between  two  solutions. 

(b)  It  is  defined  in  the  context  of  a  set  of  solutions. 

(c)  It  is  not  affected  by  the  scales  of  the  different  objectives  since  the  same 
objectives  from  different  solutions  occur  in  the  numerator  and  the 
denominator. 
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We  use  this  concept  in  the  context  of  a  set  of  non-dominated  solutions.  For  each 
solution  in  a  set  of  n  non-dominated  solutions,  it  is  proposed  to  define  the  Level  of 
Dominance  (LoD)  pt  for  a  solution  i  as  follows: 

A  =Z/L  (3-5) 

k= 1 

Thus  pj  gives  the  sum  total  of  the  degrees  of  domination  for  a  given  solution  and 
thus  a  solution  with  the  lowest  value  of  pt  is  the  least  dominated  solution  in  the  set. 

This  would  allow  us  to  choose  a  single  solution  from  a  set  of  Pareto  solutions.  There 
are  many  techniques  in  the  literature  for  design  selection  amongst  Pareto  solutions, 
but  most  of  them  require  additional  information  about  the  designs  or  some  kind  of 
preference  articulation  about  the  different  objectives.  The  method  proposed  above 
does  not  require  any  such  information  and  can  be  used  to  choose  designs  based  solely 
on  the  relative  magnitudes  of  the  objective  function  values.  Other  similar  methods 
can  be  used  as  well  (Deb,  2001). 

This  method  is  used  in  subsequent  sections  when  it  is  required  to  choose  a  single 
solution  from  a  set  of  Pareto  solutions  obtained  from  solving  intennediate 
optimization  problems  while  selecting  new  samples. 

3.4  Proposed  Cross-  Validation  Based  Approach  for  Adaptive  DOE 

This  section  describes  the  new  proposed  approach  for  single  response  adaptive 
DOE.  The  approach  is  tenned  as  Space-Filling  Cross-Validation  Tradeoff  (SFCVT) 
method.  This  approach  was  inspired  by  the  work  of  Li  and  Azarm  (2006)  and  Li 
(2007).  The  work  reported  in  the  literature  (Jin  et  ah,  2006;  Li  and  Azarm,  2006;  Li, 
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2007)  on  adaptive  design  of  experiments  focuses  on  prediction  error  and  distance 
based  criteria.  For  instance,  the  cross  validation  error  criterion  allocates  new  points  in 
the  region(s)  with  high  nonlinearity  in  the  response  space  and  can  result  in  some  areas 
of  the  design  space  with  insufficient  points  while  clustering  in  regions  with 
nonlinearity.  The  model  used  for  prediction  of  cross-validation  errors  should  be 
reasonably  robust  and  accurate  otherwise  the  model  may  predict  false  CAMM  regions 
thereby  wasting  precious  function  calls.  On  the  other  hand  any  approach  solely  based 
on  the  distance  criterion,  such  as  the  maximin  scaled  distance  method  (Jin  et  ah, 
2002)  will  most  likely  distribute  available  points  across  the  entire  design  space,  when 
intuitively  one  would  expect  more  sample  points  to  be  allocated  in  regions  with  high 
nonlinearity  or  high  cross  validation  errors.  Thus,  the  DOE  approach  should 
continuously  make  a  tradeoff  between  a  cross-validation  and  a  space  filling  (or 
distance  based)  criterion  to  sample  a  next  point.  This  leads  us  to  devise  a  two 
objective  optimization  problem  where  the  objectives  (cross-validation  and  space 
filling)  are  conflicting  since  they  compete  for  the  same  resource  (i.e.,  available 
functions  calls).  The  presumption  here  is  that  the  end-user  of  the  metamodel  is 
interested  in  the  prediction  capability  of  the  metamodel  over  the  entire  design  space 
and  not  just  local  or  global  extreme  (minimum  and  maximum)  points. 

3.4.1  Overview  of  Proposed  Approach 

We  present  two  different  problem  formulations  based  on  the  space-filling  vs. 
cross  validation  tradeoff.  The  first  formulation  termed  as  the  SO  fonnulation,  is  a 
single-objective  formulation,  in  which  the  space-filling  criterion  is  treated  as  a 
constraint  and  the  optimization  objective  is  to  find  a  point  in  the  input  space  with  the 
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highest  estimated  value  of  the  LOO  error.  In  the  second  formulation,  termed  as  MO 
formulation,  the  space-filling  criterion  is  treated  as  an  objective  along  with  cross 
validation  error.  Results  are  presented  for  both  formulations,  but  only  the  SO 
formulation  is  used  when  comparing  the  SFCVT  method  with  other  DOE  methods 
from  literature. 

The  approach  starts  out  with  an  initial  design.  Based  on  this  initial  design,  the 
LOO  errors  are  calculated  for  each  point  in  the  design  and  a  Kriging  metamodel  is 
developed  for  the  LOO  errors.  The  metamodel  is  then  used  to  predict  LOO  errors  for 
unobserved  points.  An  additional  criterion  tenned  as  “Space-filling  metric”  is 
introduced  as  well.  The  space  filling  criterion  insures  judicious  allocation  of  samples 
across  the  entire  design  space  with  emphasis  on  CAMM  regions.  The  SO  or  the  MO 
formulation  is  then  solved  to  choose  the  next  experiment.  A  step-by-step  description 
of  the  approach  along  with  a  flow  chart  is  given  in  Section  3.4.8. 

3.4.2  Choice  of  Initial  Design 

The  proposed  approach  starts  with  an  initial  design.  The  purpose  of  having  an 
initial  design  is  to  get  a  good  representation  of  the  CAMM  regions  in  the  response 
space.  The  initial  design  is  chosen  to  be  a  space-filling  design  with  a  specified 
number  of  points.  An  implicit  assumption  in  our  approach  is  that  the  initial  design 
provides  a  reasonable  representation  of  the  response  space.  Here,  a  space  filling 
design  approach  such  as  the  maximum  entropy  design  is  used;  however  any  other 
space  filling  design  method  such  as  those  mentioned  in  Chapter  2  can  be  used  as  well. 

An  important  factor  influencing  the  success  of  an  adaptive  DOE  approach  is  the 
number  of  points  in  the  initial  design.  The  number  of  points  in  the  initial  design 
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should  be  a  function  of  the  problem  dimension  and  as  such  is  very  much  problem 
dependent.  It  is  proposed  that  the  number  of  points  in  the  initial  design  should  be  half 
of  the  total  number  of  available  function  calls.  This  proposition  was  used  as  a  starting 
point  for  initial  studies  and  will  be  investigated  further.  There  is  a  possibility  that 
some  problems  may  require  fewer  number  of  points  in  the  initial  design  (and  thus 
more  adaptively  sampled  points)  and  vice-versa. 


3.4.3  LOO  Error  Prediction 

The  relative  LOO  cross-validation  error  is  used  in  the  proposed  approach  to 
explore  the  design  space  and  sample  new  points.  The  LOO  error,  by  definition,  can  be 
easily  calculated  for  the  current  points  in  the  design.  But  there  is  no  way  to  calculate 
it  for  points  not  in  the  current  design,  i.e.,  unobserved  points.  It  is  proposed  to 
develop  a  metamodel  for  the  LOO  errors  based  on  the  current  design.  This  metamodel 
can  then  be  used  to  predict  the  LOO  error  for  other  unobserved  points. 

The  metamodeling  method  used  here  is  Kriging.  As  described  in  Chapter  2,  using 
a  Gaussian  correlation,  the  LOO  error  (ewo)  can  be  correlated  in  space  using  the 
following  equation: 


eLOO(x o)  =  iu(x)  +  r  R  (y~fM) 


where,  y  = 


ewo(xi) 


eLOo(Xn ). 


(3.6) 


54 


where  ewo(x0 )  is  the  predicted  LOO  error  for  an  unobserved  points  x0,  ju(x)  is  the 

mean  of  LOO  errors  for  points  in  the  current  design,  rT  is  the  correlation  vector 
explained  in  Chapter  2  and  R  is  the  correlation  matrix  with  each  term  of  R  given  by: 

R(i,  j)  =  exp(-£  6>  |  xu  -  xfi  \pl)  (3.7) 

/= i 

In  Eq.  (3.7),  d  is  the  dimension  of  the  input  space,  and  0t  is  the  correlation 

parameter  along  each  dimension  of  the  input  space.  Based  on  studies  in  the  literature 
(Jones,  2001)  a  value  of  2  is  chosen  for  p. 

3.4.4  Space-filling  Criterion 

As  mentioned  earlier,  a  purely  cross-validation  based  approach  can  result  in 
clustered  samples  in  the  design  space.  In  order  to  avoid  this  clustering  the  use  of  a 
space  filling  criterion  is  proposed.  This  space  filling  metric  can  be  based  on  Euclidean 
distance  or  the  scaled  distance  (Jin  et  ah,  2002)  in  the  design  space. 

The  first  space  filling  metric  is  the  maximin  distance  in  the  design  space.  For  each 
experiment  in  an  existing  design  D,  we  compute  the  minimum  non-zero  distance  of 
this  point  from  all  other  points  in  D.  Thus  for  n  points,  we  have  n  non-zero  distances. 
We  compute  the  maximum  of  these  distances  and  set  the  space  filling  metric  to  be 
equal  to  one-half  of  this  maximum  value.  This  will  ensure  that  the  new  sample  point 
will  not  be  placed  close  to  existing  points  in  the  design.  This  space  filling  metric  S  is 
independent  of  the  metamodeling  technique  used.  Mathematically,  S  can  be 
represented  as: 


55 


ds{xt)  =  min(||  xi  -x .  ||2  )  ,  V  jc.  e  D  a  (/'  ^  j) 
S  =  0.5  x  max(fi?5(x;)),  Vx,.  e  D 


(3.8) 


Earlier  in  this  dissertation  we  proposed  to  use  Kriging  as  the  metamodeling 
technique. 

The  space  filling  criterion  from  Eq.  (3.8)  is  used  in  the  SO  formulation.  An 
alternative  space  filling  criterion  can  also  be  used  which  is  dependent  upon  the  use  of 
Kriging  as  the  metamodeling  technique.  This  alternative  space  filling  criterion  is 
simply  the  mean  square  prediction  error  obtained  using  the  Kriging  predictor.  The  SO 
and  the  MO  formulations  for  choosing  the  new  sample  are  described  next. 


3.4.5  Choosing  Next  Sample 

Once  a  model  for  ewo  and  the  space-filling  metric  have  been  selected,  the  next 
sample  is  chosen  such  that  it  has  the  maximum  value  of  ewo  and  is  sufficiently  far 
away  from  all  existing  points  in  the  current  design.  Mathematically  this  can  be 
represented  as: 


*B+i  =  argmaxei00 

X 

s.t.  ||jc  — jca||2  >  S  V  xk  e  D 


(3.9) 


The  problem  in  Eq.  (3.9)  tries  to  find  out  the  next  sample  as  the  one  with  the 
maximum  value  of  the  predicted  LOO  error,  eLOO ,  and  which  is  sufficiently  away 

from  all  points  in  the  current  design.  S  is  defined  in  Eq.  (3.8).  Due  to  implementation 
issues,  if  the  problem  depicted  in  Eq.  (3.9)  fails  and  cannot  be  solved,  then  the  next 
point  is  chosen  using  the  MSE  method.  The  optimization  problem  may  fail  because 
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there  are  no  feasible  solutions  remaining  implying  that  the  design  space  is  saturated. 
The  problem  may  also  fail  to  solve  due  to  implementation  issues  such  as  numerical 
rounding  and  ill-conditioned  matrices.  Note  that  all  the  candidate  points  in  the  design 
space  that  satisfy  the  constraint  are  treated  equally  with  respect  to  eLoo ■  There  is  no 
spatial  correlation  taken  into  account. 

An  alternative  formulation  using  the  mean  squared  error  as  the  space-filling 
criterion  is  proposed  as  follows: 

xM  =  arg  max  (eLOO  ,s2  )  (3.10) 

X 

9 

This  formulation  uses  the  MSE  (s  )  as  a  second  optimization  objective.  MSE 
allows  us  to  account  for  the  spatial  correlation  amongst  different  designs  in  different 
dimensions.  Moreover,  the  use  of  s~  can  be  thought  of  as  an  adaptive  space-filling 
design.  The  problem  in  Eq.  (3.10)  tries  to  find  the  next  point  such  that  it  has  the 
maximum  value  of  the  predicted  LOO  error  and  also  has  the  maximum  MSE. 

The  fonnulation  presented  in  Eq.  (3.10)  is  an  unconstrained  multiobjective 
optimization  problem  and  as  such  will  yield  a  set  of  Pareto  solutions  because  there  is 
an  inherent  tradeoff  between  the  LOO  errors  and  MSE  criteria.  All  the  Pareto  points 
can  be  chosen  as  the  next  set  of  samples  (instead  of  one  point  at  a  time).  But  this  may 
not  result  in  the  best  design.  The  challenge  is  to  select  one  single  design  at  each  stage. 
A  single  solution  can  be  chosen  based  on  one  of  the  following  procedures: 

(a)  Amongst  the  Pareto  solutions,  choose  the  one  with  maximum  value  of  ewo 
or  s  ,  i.e.,  one  of  the  extreme  Pareto  solutions.  The  decision  can  be  made 
based  on  the  relative  magnitudes  of  each  quantity  in  the  set. 
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(b)  Choose  the  solution  with  the  minimum  value  of  level  of  dominance,  which 
is  described  in  Section  3.3.3  or  first  k  solutions  with  lowest  values. 

(c)  Choose  all  the  Pareto  solutions. 

Options  (b)  and  (c)  above  are  useful  in  the  context  of  distributed  computing  when 
more  than  one  simulation  can  be  invoked  simultaneously  thereby  reducing  the  overall 
execution  time.  Choosing  one  point  per  stage  should  be  most  effective,  but  it  can  be 
very  slow  since  only  one  simulation  is  invoked  per  iteration. 

An  example  of  the  trade  off  between  the  two  objective  functions,  i.e.,  the  eLoo  and 
the  MSE  is  shown  in  Figure  3.6  for  Himmelblau’s  function  (Floudas  et  al.,  1999). 
Figure  3.6a  shows  the  actual  response  for  the  function,  3.6b  shows  the  initial  design. 
Figures  3.6c  and  3.6d  show  the  estimated  eLoo  and  the  MSE  surface.  Figure  3.6e 
shows  the  Pareto  set  obtained  by  maximizing  the  MSE  and  the  FOO  errors  and  Figure 
3.6f  shows  the  new  designs  corresponding  to  the  MSE  vs.  FOO  error  tradeoff.  It  can 
be  clearly  seen  that  there  is  a  certain  amount  of  tradeoff  between  the  two  objectives. 

In  both  cases,  the  space  filling  metric  ensures  that  the  new  samples  are 
sufficiently  away  from  the  existing  design.  Note  that  the  MSE  is  theoretically  zero  for 
points  in  the  current  design.  Thus  this  will  prevent  existing  points  from  being  sampled 
again. 

The  above  optimization  problems  can  be  solved  using  a  direct  or  exhaustive 
search  technique  or  using  a  genetic  algorithm  and  need  not  be  solved  exactly.  In  this 
dissertation,  the  exhaustive  search  was  used  for  problems  with  two  dimensions  while 
a  MOGA  (Deb,  2001)  was  used  for  problems  with  more  than  two  dimensions. 
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MSE 


Figure  3.6  (a)  Response  for  Himmelblau’s  function,  (b)  initial  design,  (c)  LOO  error 
surface,  (d)  MSE  response  surface ,  (e)  MSE  vs.  LOO  Pareto  and  (f)  new  designs  based 

on  MSE  vs.  LOO  Pareto 
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The  end-result  of  adaptive  sampling  should  result  in  a  acceptable  metamodel. 
Figures  3.7,  3.8  and  3.9  show  the  effect  of  using  the  above  two  formulations  on  the 
performance  of  the  resulting  metamodels.  In  this  example  we  choose  the  test  problem 
P3.  The  single  objective  (SO)  SFCVT  formulation  is  compared  against  the 
multiobjective  formulation.  For  the  multiobjective  formulation,  we  choose  samples 
based  on  the  level  of  dominance  (MO-LoD),  the  point  closest  to  the  origin  based  on 
L2  norm  (MO-L2)  and  the  two  extreme  (maximum  MSE  value,  MO-Max-MSE,  and 
maximum  LOO  error  value,  MO-Max-LOO)  points.  The  initial  number  of  samples 
for  each  case  was  20.  Additional  10  samples  were  added  using  the  SFCVT  method. 
Metamodels  were  developed  and  evaluated  using  a  random  sample  after  every  two 
additional  points.  The  metamodel  based  on  only  initial  design  was  also  evaluated, 
thus  giving  a  total  of  six  sample  steps  as  shown  in  the  figures. 


Figure  3.7  RMSE  for  problem  P3 
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Figure  3.8  MAE  for  problem  P3 


Figure  3.9  MAS  for  problem  P3 


Based  on  Figures  3. 7-3. 9,  it  is  observed  that  the  single  objective  forumulation  and 
the  multi-objective  forumuatoin  with  the  LoD  selection  provide  consistently  better 
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RMSE,  MAE  and  MAS  values.  It  is  also  seen  that  selecting  any  of  the  extreme  points 
does  not  result  in  a  good  metamodel.  While  this  behavior  is  strongly  exhibited  for  test 
problem  P3,  this  study  was  repeated  for  all  other  test  problems  (P1-P20)  and  it  was 
found  that  for  majority  of  the  test  problems  the  selection  of  the  extreme  points  did  not 
result  in  a  good  metamodel.  For  test  problem  P3,  we  also  observed  that  selecting  the 
next  sample  based  on  the  L2  distance  does  not  necessarily  yield  a  good  metamodel. 
This  is  because  the  L2  distance  calculation  is  affected  by  the  magnitudes  and  the 
distribution  of  the  resulting  Pareto  solutions  (MSE  vs.  LOO)  in  the  MO  formulation. 
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3.4.6  Design  Update 

The  simulation  is  invoked  for  the  sample  obtained  in  the  previous  step  and  then 
this  experiment  is  added  to  the  current  design  D.  In  this  way,  the  approach  can  sample 
as  many  points  as  desired  by  the  end-user. 

3.4.7  Stopping  Criteria 

In  real  world  engineering  design  problems,  resources  are  limited  and  generally  a 
limit  is  placed  on  the  available  number  of  function  calls  or  simulation  calls.  This  limit 
is  used  as  the  stopping  criterion  for  the  proposed  approach.  When  computational 
resources  are  readily  available,  but  still  limited,  an  elaborate  stopping  criterion  such 
as  the  steps  below  should  be  used: 

Step-1:  Sample  sufficient  number  of  random  points  and  evaluate  their  true 
response.  These  random  points  will  be  used  to  evaluate  the  metamodel 
performance. 

Step-2:  During  each  stage  of  the  proposed  approach  or  during  intennediate  stages 
(e.g.,  once  every  5  new  samples),  develop  a  metamodel  using  the  current  design 
and  evaluate  its  performance  against  the  random  sampled  points  in  Step- 1 . 

Step-3:  If  the  RMSE,  MAE,  RRMSE,  RMAE  and  MAS  values  obtained  after 
predicting  the  response  for  the  random  samples  using  the  metamodel  are 
acceptable,  stop  sampling  new  points. 

3.4.8  Step-by-Step  Description  of  Proposed  Approach 

Figure  3.10  shows  a  flow  chart  for  the  proposed  approach.  This  flow  chart 

assumes  that  the  stopping  criterion  used  is  the  maximum  number  of  available  function 
calls. 
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Figure  3.10  Flow  chart  for  the  proposed  single  response  adaptive  DOE  method 

Let  Kmax  be  the  maximum  number  of  available  function  calls.  The  different  steps 
are  as  follows: 

Step-1:  Generate  an  initial  design  D  using  the  MED  method  and  the  chosen 
number  of  initial  samples. 

Step-2:  Calculate  LOO  error  for  all  the  points  in  the  initial  design  and  build  a 
Kriging  based  metamodel  for  the  LOO  errors. 

Step-3:  Solve  the  optimization  problem  discussed  earlier  in  Eq.  (3.9)  or  (3.10)  to 
obtain  the  next  sample  point.  In  this  dissertation,  both  formulations  are  used. 

Step-4:  If  the  optimization  fails  or  if  there  are  no  feasible  solutions,  then  go  to 
Step-7;  else  go  to  Step-5. 
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Step-5:  Evaluate  the  true  response  for  the  new  point  and  add  it  to  D. 

Step-6:  If  the  number  of  points  in  D  is  equal  to  the  available  number  of  function 
calls,  then  Stop,  otherwise  continue  to  Step-2. 

Step-7:  Compute  the  next  sample  using  either  the  MED  or  the  MSE  method  and 
proceed  to  Step-5. 

The  computational  complexity  of  this  approach  is  discussed  in  Section  3.5.6  along 
with  actual  CPU  times  obtained  from  numerical  experiments. 

3.5  Application  Examples 

In  this  section,  the  proposed  approach  is  applied  to  several  numerical  test 
problems  from  a  test  suite  developed  as  a  part  of  this  dissertation.  The  problems  in  the 
test  suite  have  been  carefully  chosen  to  evaluate  the  different  aspects  of  simulation 
responses  such  as  being  highly  nonlinear  and  multimodal  in  nature  and  to 
demonstrate  the  scalability  of  the  proposed  approach.  The  details  (Sasena  et  ah,  2000; 
Jin  et  ah,  2002;  Busby  et  ah,  2007)  of  the  test  problems  are  described  in  Appendix-A. 
All  of  the  test  problems  are  of  multi-input  single-output  type  with  the  number  of 
inputs  ranging  from  two  to  eight.  The  test  problems  are  designated  as  P 1  through  P24 
for  brevity.  The  prefix  ‘P’  is  omitted  in  certain  plots  due  to  space  restrictions.  The 
Kriging  implementation  available  in  the  MATLAB  (Matlab,  2007)  DACE  Toolbox 
(Lophaven  et  ah,  2002)  was  used. 

The  proposed  approach  using  the  SO  formulation  is  compared  to  other  single 
response  DOE  approaches  from  the  literature.  The  non-adaptive  approach  used  for 
comparison  is  the  MED  (Shewry  and  Wynn,  1987;  Koehler  and  Owen,  1996)  method 
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and  two  adaptive  DOE  methods  used  for  comparison  are  the  MSD  method  (Jin  et  ah, 
2002)  and  the  ACE  method  (Li  and  Azarm,  2006;  Li,  2007)  described  earlier  in 
Chapter  2.  The  comparison  is  perfonned  on  the  basis  of  metamodel  predictions  for  a 
random  test  sample  using  a  fixed  number  of  samples  (i.e.,  fixed  number  of  function 
evaluations).  Since  random  samples  are  used  for  evaluating  the  performance  of 
resulting  metamodels,  the  numerical  experiments  are  repeated  10  times  each  to  obtain 
descriptive  statistics  for  the  results.  The  repeat  count  of  10  was  chosen  based  on 
available  computational  resources.  The  descriptive  statistics  are  presented  in  the  fonn 
of  a  box  plot  (McGill  et  ah,  1978),  also  known  as  a  box  and  whisker  plot. 

The  test  matrix  used  for  conducting  numerical  experiments  is  shown  in  Table  3.1. 
The  table  shows  the  test  problem  numbers,  the  problem  size  i.e.,  input  dimension  or 
number  of  inputs,  the  number  of  points  in  the  initial  design,  the  number  of  new 
sampled  points,  number  of  random  test  points  used  and  the  number  of  runs  conducted 
for  repeatability.  The  initial  design  was  generated  using  the  MED  method.  The  SO 
formulation  from  Eq.  (3.9)  and  the  MO  formulation  Eq.  (3.10)  were  used  to  obtain 
new  sample  points  and  the  perfonnance  of  resulting  metamodels  were  compared. 
When  comparing  the  proposed  approach  with  those  from  the  literature,  only  the  SO 
formulation  was  used.  The  metamodel  acceptability  threshold  £’MAS  was  set  to  0.1 
(i.e.,  10%)  for  the  lower  as  well  as  upper  bound. 
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Table  3.1  Test  Matrix  for  Comparison 


Test 

Problems 

Problem 

Input 

Dimension 

Initial 

Design 

New 

Samples 

Random 

Test 

Samples 

Repeat 

Runs 

P1-P16 

2 

20 

20 

2000 

10 

P17-P20 

4 

40 

40 

4000 

10 

P21-P24 

8 

80 

80 

8000 

10 

A  Factor  of  Improvement  metric  is  introduced  to  compare  metamodel 
performance  and  to  visualize  the  comparison  results  as  described  next. 

3.5.1  Metamodel  Comparison:  Factor  of  Improvement 

The  accuracy  of  two  metamodels  can  be  compared  by  calculating  the  relative 

reduction  in  RMSE  or  MAE  between  the  two  models  as  follows:  The  Factor  of 
Improvement  (Fol)  for  metamodel  RMSE  of  the  proposed  approach  over  another 
baseline  approach  is  defined  as: 

RMSE  FoI=(RMSEB„„lin,  -  RMSEmpo„d )  x  1 00/RMSEB.„,„  (3.11) 

Thus,  a  positive  Fol  indicates  that  the  ‘Proposed’  metamodel  is  superior  to  the 
‘Baseline’  metamodel.  Similarly,  Fol  for  MAE,  RRMSE  and  RMAE  can  be  defined. 
Clarke  et  al.  (2005)  propose  a  similar  metric  based  on  using  the  standard  deviation  of 
the  RMSE  instead  of  the  actual  values. 

3.5.2  Results  for  the  Proposed  Approach 

The  MAS  of  the  proposed  approach  using  the  SO  formulation  are  shown  for  the 

different  test  problems  in  Figures  3.11  and  3.12.  Figure  3.11  shows  the  MAS  values 
for  all  the  test  problems  for  one  single  numerical  experiment.  Figure  3.12  shows  the 
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descriptive  statistics  for  MAS,  obtained  after  conducting  the  numerical  experiments 
10  times. 

We  observe  from  Figure  3.11  that  16  out  of  24  test  problems  have  a  score  of  less 
than  50.  One  reason  for  this  is  that  a  fixed  correlation  was  used  for  Kriging 
metamodels.  The  type  of  test  problem  is  not  known  a  priori  and  thus  there  is  no  way 
to  choose  a  best  suited  correlation  except  by  trial  and  error. 

The  plot  shown  in  Figure  3.12  needs  some  explanation.  The  plot  is  termed  as  a 
box  plot  and  shows  the  following  five  parameters  in  order:  minimum  value,  1st 
quartile,  median,  3ld  quartile  and  the  maximum  value.  Thus  from  Figure  3.12,  the 
vertical  line  in  the  box  is  the  median  value.  Outliers  are  represented  by  the  “+”  signs 
on  either  end  of  the  minimum  or  the  maximum  values. 
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Figure  3.11  MAS  for  designs  from  the  proposed  SFCVT  approach 

Figure  3.13  shows  the  comparison  of  the  MAS  values  obtained  using  the  two 
different  formulations,  i.e.,  the  SO  and  the  MO  formulations.  The  level  of  dominance 
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criteria  was  used  to  select  one  point  at  a  time  for  the  MO  formulation.  It  can  be  seen 
that  the  two  formulations  yield  comparable  results  (MAS  within  ±5)  for  1 7  out  of  24 
test  problems.  Figure  3.14  shows  MAS  comparison  for  the  two  formulations,  when 
extreme  points  are  selected  at  each  stage.  Again,  17  out  of  24  test  problems  show 
comparable  results.  The  corresponding  RRMSE  and  RMAE  values  (not  shown  here) 
for  the  two  formulation  are  also  comparable.  The  MO  formulation  can  be  used  when 
multiple  sample  points  are  desired  per  stage,  which  can  then  be  evaluated 
simultaneously  taking  advantage  of  parallel  computing. 


P24 

- 

1 

I 

m  *  ■  ^ 

1  1 

- 

P23 

- 

1 

- 

P22 

4-  Ob 

- 

P21 

- 

♦  id 

- 

P20 

- 

(D-1 

- 

P19 

- 

i —  i — i 

- 

P18 

- 

+ 

- 

P17 

- 

i  i — i 

- 

P16 

- 

<3 

- 

P15 

Ijgp 

- 

P14 

- 

■-*U* 

- 

P13 

|if.i 

- 

P12 

,§j! 

- 

Pit 

- 

M1 

- 

P10 

- 

- 

P9 

■1 

- 

P8 

- 

P7 

- 

- 

P6 

- 

f 

- 

P5 

- 

# 

- 

P4 

- 

- 

P3 

- 

. 

- 

P2 

fife# 

- 

PI 

- 

i 

1 

1 _  1  I  1 

1 

T  + 

- 

20 

30 

40  50  60  70 

80 

90  100 

MAS  for  SFCVT 

Figure  3.12  MAS  for  the  proposed  SFCVT  approach 
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Figure  3.13  MAS  of  the  SO  and  MO  formulations,  using  LoD  as  selection  criterion  in 

MO  formulation 
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The  RMSE,  MAE,  RRMSE  and  RMAE  plots  results  for  the  proposed  approach 
are  shown  in  the  following  sections  where  the  SFCVT  method  is  compared  with  three 
methods  from  the  literature.  The  SO  formulation  was  used  for  this  comparison.  In 
cases  where  results  are  presented  in  the  fonn  of  box  plots,  conclusions  are  drawn 
based  on  the  median  values  shown  in  the  box  plot.  In  some  cases,  when  all  values  in  a 
box  plot  are  very  close,  the  entire  plot  appears  as  only  one  horizontal  or  vertical 
(based  on  plot  orientation)  line. 

3.5.3  Comparison  with  MED  Method 

In  this  section,  the  proposed  approach  is  compared  with  the  MED  approach 
(Shewry  and  Wynn,  1987).  In  order  to  ensure  a  fair  comparison,  the  number  of  initial 
designs  and  the  total  number  of  designs  was  fixed.  Numerical  experiments  were 
carried  out  as  per  the  test  matrix  presented  in  Table  3.1.  The  results  for  MAS  for  the 
MED  method  and  the  proposed  SFCVT  method  are  shown  in  Figure  3.15.  In  this 
comparison  it  can  be  observed  that  the  proposed  approach  outperforms  (i.e.,  yields 
higher  MAS)  the  MED  method  in  terms  of  MAS  in  1 1  of  24  test  problems.  Figure 
3.13  shows  the  box  plot  for  MAS  values  for  the  MED  method.  The  horizontal  length 
of  the  boxes  is  small,  indicating  a  consistent  performance  of  the  MED  method  on  the 
random  sample  points  used  for  verification.  Upon  further  investigation  it  was  found 
that  the  default  correlation  function  used  in  the  SFCVT  method  for  the  LOO  error 
prediction  was  not  suitable  i.e.,  it  did  not  give  a  good  prediction  of  the  LOO  errors. 
Changing  the  correlation  to  an  exponential  function  (see  Chapter  2)  improved  the 
results  dramatically  for  test  problem  P4,  P5,  P9  and  P14-P16. 


71 


Figures  3.17,  3.18,  3.18  and  3.20  show  the  RMSE  Fol,  MAE  Fol,  RRMSE  Fol 
and  RMAE  Fol  of  the  proposed  SFCVT  method  over  the  MED  method.  The  proposed 
approach  was  better  in  terms  of  RMSE  for  15,  in  terms  of  MAE  for  16,  in  terms  of 
RRMSE  for  16  and  in  terms  of  RMAE  for  18  out  of  24  test  problems.  These  numbers 
were  obtained  based  on  the  median  values  from  Figures  3.17-3.20.  We  can  define  the 
tenn  “strictly  better”  for  Fol,  when  all  the  values  represented  by  the  box  plot  are 
greater  than  0.  Thus,  the  proposed  approach  performed  strictly  better  than  the  MED 
method  in  terms  of  RMSE  for  10,  in  terms  of  MAE  for  12,  in  terms  of  RRMSE  for  1 1 
and  in  terms  of  RMAE  for  13  out  of  24  test  problems.  The  proposed  approach  faired 
the  worst  on  problem  P10.  Upon  further  investigation,  it  was  found  that  the 
metamodel  for  the  LOO  errors  was  not  sufficiently  accurate  and  as  a  result  only  the 
space-filling  methods  (i.e.,  MED  and  MSD)  showed  an  acceptable  metamodel. 
Moreover,  it  was  observed  that  for  problem  P10,  the  LOO  error  metamodel  was 
sensitive  to  the  initial  maximum  entropy  design  sample  and  as  such  never  yielded  a 
good  model  for  the  LOO  errors  in  intermediate  stages,  ultimately  resulting  in  an 
inferior  design.  Further  studies  on  initial  design  were  conducted  for  problem  P10.  It 
was  found  that  if  LHC  design  is  used  as  the  initial  design  or  even  if  the  prior  variance 
and  correlation  parameters  were  changed  for  the  MED  method  to  generate  initial 
design,  the  resulting  metamodel  has  a  better  performance  than  that  obtained  using  the 
MED  method. 
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Figure  3.15  MAS  comparison  for  MED  vs.  SFCVT 
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Figure  3.16  MAS  for  MED  approach 
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Figure  3.17  RMSE  Fol  of  proposed  approach  over  MED  method 
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Figure  3.18  MAE  Fol  of  proposed  approach  over  MED  method 
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RRMSE  Fol  over  MED 


Figure  3.19  RRMSE  Fol  of  proposed  approach  over  MED  method 


Figure  3.20  RMAE  Fol  of  proposed  approach  over  MED  method 
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In  Figures  3.17-3.20,  the  width  of  the  box  plots  for  some  of  the  test  problems  is 
very  long,  e.g.,  P21  in  Figure  3.17.  This  indicates  that  the  number  of  repeat  runs  was 
not  sufficient  and  more  numerical  experiments  need  to  be  conducted. 

The  problem  P10  is  further  analyzed  here.  Figure  3.21a  and  3.21b  show  the 
design  generated  using  the  MED  and  the  SFCVT  methods,  respectively.  The 
corresponding  verification  error  (i.e.,  absolute  value  of  the  difference  between  actual 
and  predicted  response)  plots  are  shown  in  Figures  3.21c  and  3.2 Id  for  MED  and 
SFCVT  respectively.  For  this  case,  the  MAS  for  SFCVT  is  much  lower  than  that  for 
MED.  We  now  tweak  the  correlation  parameter  6  used  in  the  MED  method.  The  new 
designs  are  shown  in  Figure  3.22a  and  3.22b.  The  corresponding  error  scatters  are 
shown  in  Figures  3.22c  and  3.22d  for  MED  and  SFCVT  respectively.  Note  that  even 
though  the  overall  MED  design  hasn’t  changed  significantly,  the  design  generated 
using  the  SFCVT  method  is  drastically  different  and  its  performance  is  now 
comparable  to  that  of  the  MED  method.  This  shows  that  for  certain  problems,  the 
SFCVT  approach  is  sensitive  to  the  initial  design. 
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Figure  3.21  Problem  P10,  (a)  MED  design,  (b)  SFCVT  design,  (c)  absolute  error  using 
MED  and  (d)  absolute  errors  using  SFCVT 
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Figure  3.22  Problem  P10  with  modified  initial  design,  (a)  MED  design,  (b)  SFCVT 
design,  (c)  absolute  error  using  MED  and  (d)  absolute  errors  using  SFCVT 


3.5.4  Comparison  with  MSD  Method 

In  this  section,  the  proposed  approach  is  compared  with  the  MSD  method  (Jin  et 
ah,  2002).  Again,  the  number  of  initial  designs  and  the  total  number  of  designs  was 
fixed.  The  results  are  shown  in  Figures  3.23-3.28. 

From  Figure  3.23,  we  see  that  the  proposed  approach  outperforms  the  MSD 
method  in  only  6  out  of  24  test  problems.  The  performance  improvement  of  the  MSD 
method  over  the  proposed  SFCVT  method  in  terms  of  MAS  is  small,  except  for 
problems  P7,  P9  and  P10.  Figure  3.23  shows  the  box  plots  for  MAS  values  obtained 
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from  the  MSD  method.  Figures  3.25,  3.26,  3.27  and  3.28  show  the  RMSE  Fol,  MAE 
Fol,  RRMSE  Fol  and  RMAE  Fol  of  the  SFCVT  approach  over  MSD  method.  Based 
on  the  median  values,  the  SFCVT  method  performs  better  than  MSD  method  in  12 
out  of  24  problems  for  RMSE,  1 1  out  of  24  for  MAE,  13  out  of  24  for  RRMSE  and 
13  out  of  24  for  RMAE. 


The  LOO  error  model  issue  for  problem  P10,  described  in  previous  section  is 
present  in  this  case  as  well.  Only  the  MSD  and  the  MED  methods  work  well  on 
problem  P10. 
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Figure  3.23  MAS  comparison  for  MSD  vs.  SFCVT 
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Figure  3.24  MAS  for  MSD  method 
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Figure  3.25  RMSE  Fol  of  proposed  approach  over  MSD  method 
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Figure  3.26  MAE  Fol  of  proposed  approach  over  MSD  method 
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Figure  3.27  RRMSE  Fol  of  proposed  approach  over  MSD  method 
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Figure  3.28  RMAE  Fol  of  proposed  approach  over  MSD  method 


3.5.5  Comparison  with  ACE  Method 

The  proposed  approach  is  also  compared  with  the  ACE  method  (Li  and  Azarm, 
2006;  Li,  2007).  Due  to  programming  issues  with  the  available  ACE  implementation, 
the  comparison  could  only  be  carried  out  for  problems  with  two  input  variables. 
Ligure  3.29  shows  the  comparison  of  MAS  for  the  SFCVT  and  the  ACE  method.  It 
can  be  observed  from  Figure  3.29  that  the  proposed  SFCVT  approach  has  better  MAS 
than  that  of  ACE  method  for  11  out  of  16  test  problems.  The  other  Fol  plots  are 
shown  in  Figures  3.30  -  3.33.  Again,  based  on  median  values,  the  SFCVT  method 
performs  better  than  the  ACE  method  in  1 1  out  of  16  problems  in  terms  of  RMSE,  1 1 
in  terms  of  MAE,  10  in  terms  of  RRMSE  and  1 1  in  terms  of  RMAE. 
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Figure  3.29  MAS  comparison  for  ACE  vs.  SFCVT 
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Figure  3.30  RMSE  Fol  of  proposed  approach  over  ACE  method 


83 


P16 

- 

P15 

- 

P14 

- 

P13 

- 

P12 

- 

Pll 

- 

a  P10 

- 

bier 

hs 

SO 

- 

o 

tPr 

h d 
00 

Tes 

hs 

- 

P6 

- 

P5 

- 

P4 

- 

P3 

- 

P2 

- 

PI 

-150 

Tfpf 

CH 


m 


i— £j 


d — 1 

HD 

i§ 


_ i _ 

-100 


_ I _ L_ 

-50  0 

MAE  Fol  over  ACE 


50 


I 


Figure  3.31  MAE  Fol  of  proposed  approach  over  ACE  method 
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Figure  3.32  RRMSE  Fol  of  proposed  approach  over  ACE  method 
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The  proposed  approach  is  also  applied  to  a  real-world  engineering  design  problem 
which  focuses  on  the  research  and  development  of  the  next  generation  of  air-cooled 
heat  exchangers.  The  details  and  the  results  for  this  engineering  example  are 
discussed  in  Chapter  5. 

3.5.6  Effect  of  Initial  Design 

In  this  section  we  briefly  study  the  effect  of  the  number  of  points  in  the  initial 
design.  For  any  given  problem,  there  are  an  optimal  number  of  points  that  should  be 
present  in  the  initial  design.  While  this  number  cannot  be  arrived  apriori  for  black  box 
functions,  in  this  section  we  study  the  effect  of  the  number  of  points  in  the  initial 
design  for  a  fixed  number  of  total  available  function  calls.  The  effect  of  the  choice  of 
the  initial  design  can  also  be  studied. 
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The  MED  method  is  used  to  generate  initial  samples.  The  total  number  of  sample 
points  is  fixed  at  (20 xd),  where  d  is  the  problem  dimension.  Numerical  experiments 
are  conducted  starting  with  5 d  and  10 d  as  the  number  of  samples  in  the  initial  design. 
The  verification  of  the  resulting  metamodels,  expressed  using  MAS,  RRMSE  and 
RMAE  are  shown  in  Figures  3.34-3.36  below. 
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Figure  3.34  MAS  comparison  for  5 d  and  10 d  points  in  initial  design 


Figure  3.34  shows  that  for  certain  problems  (e.g.,  problems  5,7  and  10)  the  fewer 
number  of  points  in  the  initial  design  results  in  a  better  metamodel,  thus  pointing  out 
that  more  adaptive  samples  are  necessary  to  capture  function  behavior  in  the  CAMM 
regions.  The  metamodels  for  other  problems  (e.g.,  4,  5,  9  and  17)  perform  better  with 
a  higher  number  of  samples  in  the  initial  design.  This  is  because  a  lower  number  of 
initial  samples  are  not  sufficient  to  capture  all  the  CAMM  regions  in  the  response 
space.  Similar  behavior  is  observed  for  RRMSE  and  RMAE  as  shown  in  Figures  3.35 
and  3.36  respectively. 
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Figure  3.35  Metamodel  RRMSE  for  5 d  and  1  Qd  points  in  initial  design 


Figure  3.36  Metamodel  RMAE  for  5 d  and  10 d  points  in  initial  design 
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3.5.7  Scalability  to  Higher  Dimensions 

The  SFCVT  method  was  also  applied  to  test  problems  with  higher  dimensions. 
Four  test  problems  with  scalable  closed  form  expressions  were  chosen.  The  single 
objective  formulation  was  used.  The  16  and  24  dimensional  versions  of  test  problems 
P3  (Dixon-Price  function)  and  P5  (Griewank  function)  were  used.  The  test  problems 
are  described  in  Appendix-A.  The  MAS  values  were  compared  with  those  obtained 
using  the  MED  method  and  the  MSD  method.  The  comparison  is  shown  in  Figure 
3.37.  The  MAS  values  are  very  comparable  for  test  problem  P3  for  16  dimensions, 
but  for  other  test  problems,  the  SFCVT  method  performs  worse  than  the  MED  and 
the  MSD  method.  This  can  be  explained  by  the  fact  that  the  test  problems  are  highly 
multi-modal  as  such  the  method  gets  trapped  into  a  subset  of  the  CAMM  regions.  For 
such  functions  (e.g.,  Griewank  function),  a  space-fdling  approach  is  better  suited.  It 
should  be  pointed  out  that  the  function  behavior  is  assumed  to  be  unknown. 


70 


P3,d=16  P5,  d=16  P3,  d=24  P5,  d=24 

Test  Problem 


Figure  3.37  MAS  comparison  when  SFCVT  method  is  applied  to  higher  dimension 
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3.5.8  Comments  on  Computational  Effort 

The  computation  time  for  the  proposed  SFCVT  approach  is  dominated  by  the 

LOO  error  calculation  and  the  optimization  problem  which  is  solved  to  arrive  at  the 
next  sample.  Note  that  in  the  LOO  error  calculation,  the  correlation  parameters  are  re- 
estimated  for  every  point  in  the  current  design. 

In  general,  for  the  same  number  of  function  calls,  any  adaptive  DOE  approach 
will  require  more  computation  time  than  a  non-adaptive  approach.  In  order  to 
generate  the  same  number  of  new  samples,  the  SFCVT  method  requires  on  average 
twice  as  much  time  as  the  MED,  MSD  or  the  ACE  method.  While,  for  the  test  case 
considered,  the  comparison  is  made  for  the  same  number  of  function  calls,  but  at  the 
end  of  the  sampling  process  the  resulting  metamodels  have  differing  accuracies.  For 
example,  it  is  observed  that  the  SFCVT  method  performs  better  (in  terms  of  its 
accuracy  for  the  same  number  of  function  calls)  on  test  problems  with  non-symmetric 
CAMM  regions  (e.g.,  test  problems  PI,  P3,  P6,  P8),  where  as  for  test  problems  with 
multiple  distributed  and  symmetric  CAMM  regions  (e.g.,  test  problems  P4,  P7,  P10)  a 
space  fdling  (i.e.,  non-adaptive  approach)  is  better  suited.  As  such,  for  test  problems 
with  non-symmetric  CAMM  regions,  SFCVT  method  will  require  less  computational 
effort  to  arrive  at  a  globally  accurate  metamodel. 
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3.6  Empirical  Procedure  for  Function  Call  Allocation 


In  this  section  we  propose  an  empirical  procedure  for  resource  allocation  for 
single  response  adaptive  DOE  method.  The  resource  refers  to  the  available  number  of 
function  calls. 

3.6.1  Problem  Statement 

In  practical  applications,  the  number  of  available  function  calls  is  limited.  These 
function  calls  need  to  be  used  for  DOE  phase  and  metamodel  verification.  Thus  a 
systematic  method  is  required  for  the  allocation  of  function  calls,  either  for 
metamodel  development  (i.e.,  adaptive  sampling)  or  for  random  verification. 

This  problem  can  be  greatly  simplified  if  some  information  about  the  behavior  of 
the  function  is  known  before  hand.  As  with  the  assumptions  earlier,  we  assume  that 
no  information  about  the  function  is  available  a  priori. 

3.6.2  Proposed  Empirical  Procedure 

The  proposed  empirical  procedure  is  based  on  the  notion  that  as  more  and  more 
samples  are  added  to  a  design,  the  performance  of  the  metamodel  does  not  improve 
proportionately.  In  other  words,  the  metamodel  is  already  at  its  best  (i.e.,  no 
improvement  in  MAS,  RMSE,  MAE  etc.)  and  new  samples  do  not  add  any  new 
information  to  it. 

The  procedure  starts  out  with  a  set  of  inputs  from  the  designer.  The  inputs 
required  are  the  maximum  available  function  calls  (Kmax)  ,  minimum  number  of 
points  to  be  used  for  verification  (Nv)  and  the  minimum  required  MAS  (MAScrit).  It  is 
implicit  that  the  procedure  has  access  to  the  simulation  tool  to  carry  out  the  required 
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function  evaluations.  Let  d  be  the  problem  dimension,  i.e.,  the  number  of  inputs 
variables,  D  be  the  current  design.  Note  that  D  is  initially  empty. 

The  different  steps  are  explained  below. 

Step-1 :  Accept  the  inputs  Fmax  and  Nv.  Set  counter  k  =0. 

Step-2:  Generate  random  test  points  and  evaluate  their  true  response.  Update 
counter  k. 

Step-3:  Generate  initial  design  using  MED  or  other  suitable  method  and  evaluate 
the  true  responses.  The  number  of  points  in  the  initial  design  is  set  to 

Ninit  =  2  +  ( J  + 1)  x  ( J  +  2)  /  2  (3.12) 


The  reason  for  using  Eq.  (3.12)  is  as  follows:  the  second  term  in  Eq.  (3.12)  is  the 
number  of  points  required  in  a  second  order  polynomial  regression  equation  which  is 
used  in  the  response  mean  calculation  in  Kriging.  The  first  tenn  (i.e.,  2)  ensures  that 
during  LOO  error  calculation,  the  resulting  regression  is  a  valid  least-squares  system. 

Step-4:  Use  the  adaptive  DOE  method  proposed  earlier  to  generate  d  additional 
points. 

Step-5:  Develop  a  metamodel  using  the  current  design  and  then  evaluate  its 
performance  using  the  random  points  generated  in  Step-2.  Compare  its  performance 
with  metamodel  in  previous  iteration.  If  (k+d)  >  Kmax,  then  STOP.  If  this  is  the  first 
iteration,  then  go  to  Step-9.  If  the  current  MAS  is  greater  than  the  previous  MAS  and 
if  the  current  MAS  is  greater  than  MAScrit,  then  go  to  Step-6  else  go  to  Step-7. 
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Step-6:  If  (k  +d)  >  Kmax),  then  STOP,  else  Sample  d  new  points  randomly  and  add 
them  the  to  the  verification  design.  Go  to  Step- 10. 

Step-7:  Calculate  the  LOO  errors  and  the  standardized  residual  for  the  current 
design.  If  the  number  of  points  with  a  standardized  residual  greater  than  3,  is  more 
than  two,  go  to  Step-8,  else  go  to  Step-9.  The  standardized  residual  is  just  the  ratio  of 
the  cross-validation  error  and  the  mean  error  predicted  by  Kriging. 

Step-8:  Sample  d  new  points  using  the  MED  or  other  suitable  space-filling 
method  and  compute  the  corresponding  responses.  Update  D  with  the  newly 
generated  d  points.  Go  to  Step- 10. 

Step-9:  Sample  d  new  points  using  the  adaptive  DOE  method  and  compute  their 
response.  Update  D  with  the  newly  generated  d  points. 

Step- 10:  Update  counter  k,  k  =  k  +  d.  If  (k  <  Kmax),  go  to  Step-5  else  STOP. 

3.6.3  Numerical  Experiments 

The  procedure  is  applied  to  the  test  problems  PI  through  PI 6.  The  numerical 
experiments  are  conducted  20  times,  using  the  following  parameters,  Kmax  =  60,  Nv  = 
20,  MAScrit  =  50.  The  MED  method  was  used  to  generate  the  initial  design  and  the 
adaptive  DOE  method  proposed  in  section  3.4  is  used  to  adaptive  sample  new  points. 
The  results  are  shown  in  the  form  of  a  box  plot.  Figure  3.38  shows  the  distribution  of 
the  number  of  points  used  in  developing  the  metamodel.  In  10  out  of  16  test 
problems,  all  available  points  were  used  for  adaptively  sampling  new  points.  This  is 
because,  the  criteria  for  acceptance  of  the  metamodel  was  too  strict  (MAScrit  =  50). 
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Figure  3.41  shows  the  number  of  points  used  in  verification.  The  sum  of  points  in 


Figures  3.40  and  3.41  for  each  case  must  be  equal  to  Kmax. 
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Figure  3.38  Points  used  for  developing  metamodel 
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Figure  3.39  Points  used  for  verification 
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Figure  3.40  Resulting  MAS  for  test  problems 


Figure  3.40  shows  the  MAS  for  the  resulting  metamodels.  Again,  it  is  observed 
that  the  resulting  metamodels  fall  short  of  the  prescribed  MAS  criteria  on  several  test 
problems.  As  a  result,  the  procedure  never  samples  additional  points  for  verification. 

3.7  Summary 

An  adaptive  DOE  approach  for  single  response  deterministic  simulations  is 
developed  based  on  the  tradeoff  between  cross-validation  errors  and  space-filling 
criterion.  The  proposed  approach  is  based  on  the  premise  that  at  each  stage  during  an 
adaptive  DOE  procedure,  a  tradeoff  has  to  be  made  between  allocating  new  sample  in 
a  region  of  the  design  space  that  is  sparsely  sampled  and  allocating  a  new  sample  in  a 
region  where  the  response  is  estimated  to  be  multi-modal,  non-smooth  or  noisy.  The 
proposed  approach  uses  the  concept  of  cross-validation  error,  specifically  leave-one- 
out  error  to  probe  the  response  space  for  sensitive  regions.  It  also  uses  a  space  filling 
metric  to  identify  the  sparsely  populated  regions  of  the  design  space.  Two  different 
space  filling  metrics  are  used.  First,  the  concept  of  maximin  distance  is  used;  wherein 
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all  the  dimensions  of  the  input  space  are  treated  equally  and  only  a  Euclidean  distance 
is  used  as  a  criterion.  This  gives  us  the  SO  formulation  for  the  DOE  problem.  Second, 
the  MSE  is  used  as  a  space  filling  metric  resulting  in  a  multiobjective  optimization 
problem,  termed  as  the  MO  formulation.  The  advantage  of  MSE  is  that  the  distances 
in  each  dimension  are  weighted  and  the  use  of  this  metric  is  suitable  when  Kriging  is 
used  as  the  metamodeling  approach.  The  maximin  distance  metric  can  be  used  with  or 
without  the  Kriging  metamodeling  approach.  The  MO  formulation  of  the  proposed 
approach  can  provide  one  or  more  samples  at  each  stage.  This  has  the  advantage  of 
being  able  to  sample  more  points,  hence  more  information  from  the  design  and 
response  space,  in  the  same  amount  of  clock-time.  This  is  especially  useful  when  the 
approach  is  used  along  with  a  parallel  computing  framework. 

The  proposed  approach  is  applied  to  suite  of  numerical  and  engineering  test 
problems  comprising  of  24  single  response  functions  from  the  literature.  When 
comparison  the  results  obtained  using  the  SO  and  the  MO  formulation,  it  was  found 
that  they  both  yield  comparable  results.  The  SO  formulation  is  then  compared  with 
three  methods  from  the  literature.  The  results  show  that  the  proposed  approach 
outperforms  the  maximum  entropy  design  (MED)  method  and  two  additional 
adaptive  DOE  methods,  the  MSD  and  the  ACE  method  from  the  literature,  in  most  of 
the  test  problems  with  regards  to  MAS,  RMSE  and  MAE.  Numerical  experiments  are 
conducted  several  times  to  gather  statistical  results. 

The  proposed  approach  relies  on  a  good  set  of  starting  designs  that  represent  the 
design  space  reasonably  well.  Based  on  the  numerical  experiments  it  is  observed  that 
starting  with  an  excessive  number  of  initial  designs  (thereby  less  adaptive  designs  for 
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fixed  number  of  function  calls)  actually  results  in  an  inferior  metamodel.  It  is  also 
observed  that  in  certain  cases,  the  starting  design  may  never  yield  a  good  metamodel 
for  the  LOO  errors,  causing  the  approach  to  fail. 

In  the  next  chapter  an  adaptive  DOE  method  for  multiresponse  deterministic 
simulations  is  presented. 
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Chapter  4  Multiresponse  Adaptive  Design  of  Experiments 

4.1  Introduction 

This  chapter  proposes  a  new  adaptive  DOE  method  for  multiresponse 
deterministic  computer  simulations. 

Many  real  world  engineering  simulations  are  characterized  by  multiple  inputs  and 
multiple  outputs.  The  simulation  models  are  complex,  generally  iterative  in  nature 
and  yield  all  output  values  in  a  single  execution,  i.e.,  when  a  simulation  is  executed 
with  a  set  of  inputs,  all  the  outputs  are  calculated  simultaneously.  These  outputs  or 
responses,  in  the  context  of  design  optimization,  may  be  objective  values,  constraint 
values,  design  attributes  or  simply  other  outputs  of  interest  from  the  simulation.  Naive 
application  of  approximation-assisted  optimization  to  such  a  simulation  would 
involve  developing  a  metamodel  for  each  response  separately.  In  other  words,  for 
each  output  an  adaptive  DOE  approach  is  used  to  create  a  set  of  samples.  These 
samples  may  not  be  the  same  for  each  output,  thus  wasting  a  large  number  of  function 
calls.  In  order  to  alleviate  this  problem,  multiresponse  approximation  models  (Li  et 
ah,  2006;  Li,  2007;  Romero  et  ah,  2006)  have  been  proposed  in  the  literature.  As 
opposed  to  single  response  approximation  methods,  multiresponse  approximation 
methods,  as  the  name  implies,  try  to  simultaneously  achieve  metamodels  for  multiple 
outputs  using  a  minimum  number  of  function  calls.  In  this  chapter,  the  problem  of 
adaptive  DOE  for  multiple  responses  is  posed  as  a  multiobjective  optimization 
problem,  where  there  is  a  tradeoff  between  allocating  samples  for  improving  the 
metamodel  performance  for  one  response  vs.  another.  The  proposed  approach  is 
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termed  as  Multiresponse  Space-filling  Cross-Validation  Tradeoff  (MSFCVT)  method. 
It  is  shown  that  the  proposed  approach  performs  better  than  the  previous  methods 
presented  in  the  literature  when  applied  to  the  same  test  problems. 

The  remainder  of  this  chapter  is  organized  as  follows:  Section  4.2  provides  a 
detailed  literature  review  on  existing  techniques  for  multiresponse  approximation. 
Section  4.3  describes  the  proposed  approach.  Section  4.4  demonstrates  the 
applicability  and  scalability  of  the  proposed  approach  by  applying  it  to  a  range  of  test 
problems  from  the  literature  and  comparing  its  performance  against  previous  DOE 
methods. 

4.2  Literature  Review 

There  are  two  aspects  in  multiresponse  approximation.  First  is  the  simultaneous 
DOE  for  multiple  responses.  The  DOE  method  must  sample  new  points  such  that  they 
improve  the  performance  of  the  metamodels  for  all  responses  at  the  same  time. 
Second  is  the  metamodel  development,  in  which  we  can  develop  one  metamodel  for 
each  response  or  we  can  develop  a  multiresponse  metamodel  for  all  responses 
simultaneously.  In  this  section  brief  overview  of  metamodeling  techniques  for 
multiresponse  approximation  is  presented  followed  by  a  review  of  adaptive  DOE 
techniques.  Even  though  the  DOE  is  the  first  step  in  approximation,  we  first  present 
the  metamodeling  techniques  and  then  the  DOE,  since  DOE  is  the  focus  of  this 
chapter. 

As  discussed  in  Chapter  3,  numerous  approximation  techniques  have  been 
introduced  in  the  literature  for  single  response  approximation.  The  underlying 
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assumption  of  these  metamodeling  techniques  is  that  the  engineering  design 
simulation  generates  a  single  response  i.e.,  scalar  output.  Many  complex  engineering 
simulations  yield  multiple  responses  in  a  single  execution,  i.e.,  for  a  single  set  of 
inputs  a  single  set  of  outputs  (or  responses)  is  generated.  One  can  treat  each  of  these 
responses  as  being  an  output  of  a  single  simulation  and  then  construct  a  DOE  and  a 
metamodel  for  each  of  the  responses,  which  is  a  fairly  common  approach. 

Some  research  has  been  conducted  to  handle  multiple  responses  in  a 
computationally  efficient  manner,  though  the  literature  on  multiresponse  adaptive 
design  of  experiments  is  sparse.  The  problem  of  simultaneously  handling  multiple 
responses  has  been  addressed  for  regression  methods  by  combining  the  multiple- 
responses  into  a  “value”  function  (Chiao  and  Hamada,  2001;  Derringer  and  Suich, 
1980).  An  extension  of  Kriging,  tenned  as  Cokriging  (Ver  Hoef  and  Cressie,  1993) 
has  been  proposed  to  handle  multiple  responses.  Romero  et  al.  (2006)  propose  a 
multiresponse  meta-modeling  approach  based  on  the  use  of  covariance  matrix  of  the 
responses.  This  is  very  similar  to  the  Kriging  method,  where  the  predicted  output  is 
now  a  vector  with  more  than  one  element.  They  also  propose  an  extension  to  the 
MED  method,  tenned  as  the  modified  MED  (MMED)  method,  to  handle  multiple 
responses.  The  covariance  matrix  in  the  MMED  method  is  based  on  the 
multiresponse  covariance  function  (Romero  et  al.,  2006)  of  the  underlying 
metamodel.  A  dependent  metamodeling  (DMM)  approach  (Li  et  al.,  2006)  was  also 
proposed,  wherein  instead  of  developing  a  metamodel  in  which  each  response  is  a 
function  of  the  input  design  variables,  the  response  is  assumed  to  be  a  function  of  the 
input  designs  variables  as  well  as  other  responses.  This  results  in  a  system  of 
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nonlinear  equations  (metamodels)  that  needs  to  be  solved  to  evaluate  the  multiple 
responses.  The  DMM  technique  implies  that  the  DOE  for  each  response  should  be 
constructed  such  that  it  includes  the  input  design  variables  as  well  as  the  other 
responses.  Based  on  these  inferences,  Li  et  al.  (2006)  also  include  an  extension  of  the 
maximum  entropy  design  method  to  include  the  response  space.  We  term  this 
approach  as  the  MED  for  dependent  metamodel.  Li  et  al.  (2006)  applied  their 
approach  to  examples  with  two  input  variables  and  two  responses.  Further  testing  of 
the  method  by  Li  et  al.  (2006)  is  necessary  to  evaluate  the  general  applicability  of 
their  approach. 

Romero  et  al.  (2006)  proposed  an  adaptive  DOE  approach  for  multiresponse 
approximation.  Their  approach  is  similar  to  the  Kriging  approach,  but  the  covariance 
function  is  a  multivariate  Gaussian  function  which  accounts  for  the  multiple 
responses.  The  resulting  covariance  matrix  includes  contributions  from  the  multiple 
responses  as  well.  To  complement  this  metamodeling  technique,  they  proposed  two 
adaptive  sampling  schemes.  The  first  technique  is  an  extension  of  the  maximum 
entropy  sampling  method  (i.e.,  MMED),  wherein  the  augmented  covariance  matrix  is 
replaced  by  the  new  covariance  matrix  from  Romero  et  al.’s  (2006)  method.  Second, 
the  Kleijnen  and  van  Beers  (2004)  approach  of  maximum  cross-validation  variance 
(MCVV)  is  extended  to  include  multiple  responses.  In  single  response  case,  the 
Kleijnen  and  van  Beers  (2004)  approach  finds  the  next  sample  as  the  one  that 
maximizes  the  cross-validation  variance  while  in  the  multiple  response  case,  Romero 
et  al.  (2006)  propose  to  use  the  trace  of  the  cross-validation  covariance  matrix  for  the 
multiple  responses.  The  trace  calculation  and  hence  the  approach  requires  that  the 
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responses  be  scaled  properly  in  order  to  be  effective.  Romero  et  al.  (2006)  applied 
their  approaches  to  3  test  problems  (with  2  inputs  -  1  response)  and  for  multiresponse 
case;  they  selected  two  single  response  test  problems  to  mimic  multiple  responses. 
Their  studies  found  that  the  modified  MED  method  outperformed  the  MCVV  method. 
They  used  3  test  problems  in  their  numerical  experiments  and  the  conclusion  was  that 
optimal  Latin-hypercube  (LHC)  as  initial  sampling  technique  along  with  their 
proposed  adaptive  MED  method  results  in  metamodels  that  are  more  robust  and 
accurate  or  better  than  other  sampling  approaches.  The  MCVV  method  needs  to  be 
applied  to  additional  test  problems  with  multiple  inputs  and  outputs  to  evaluate  its 
general  applicability. 

Based  on  the  literature  so  far,  the  existing  approaches  for  multiresponse  adaptive 
DOE  suffer  from  the  following  drawbacks:  (a)  scaling  and/or  normalization  of  inputs 
and  responses,  which  is  not  very  intuitive  especially  when  the  design  problem  is  not 
known  before  hand,  (b)  the  approaches  have  used  multiresponse  metamodels  and  no 
comparison  has  been  made  to  single  response  metamodels  with  multiresponse  DOE 
(c)  general  applicability  is  not  known,  i.e.,  the  approaches  have  not  been  applied  to  a 
broader  set  of  numerical  or  engineering  design  problems  with  more  than  two  inputs 
and  responses.  In  this  chapter  a  new  adaptive  DOE  method  for  multiple  responses  is 
proposed  which  extends  the  LOO  error  and  the  space  filling  criteria  defined  in 
Chapter  3  to  multiple  responses.  The  proposed  multiresponse  approach  is  applied  to  a 
total  of  10  test  problems  from  the  literature.  Two  of  these  test  problems  are  from  Li  et 
al.  (2006)  and  Romero  et  al.  (2006)  and  are  used  for  comparison  against  published 
literature. 
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4.3  Multiresponse  Adaptive  DOE  Approach 


In  this  section  we  provide  the  details  of  the  proposed  approach.  As  mentioned  in 
Chapter  3,  an  adaptive  DOE  approach  has  to  make  a  tradeoff  between  allocating 
experiments  in  sparsely  sampled  areas  versus  allocating  experiments  in  regions 
corresponding  to  the  CAMM  regions  of  the  response  space.  In  Chapter  3,  these  two 
aspects  were  handled  using  a  LOO  error  metric  and  a  space-filling  criterion.  In  the 
proposed  approach  for  this  chapter  the  same  concept  is  extended  to  multiple 
responses.  In  a  multiresponse  case,  the  tradeoff  is  not  only  between  the  LOO  error 
and  the  space-filling  criteria,  but  also  between  the  LOO  errors  for  the  individual 
responses. 

The  different  steps  in  the  proposed  approach  are  described  in  the  following 
sections.  A  detailed  flow  chart  of  the  entire  procedure  is  explained  in  Section  4.3.8. 

4.3.1  Overview  of  the  Proposed  Approach 

The  proposed  approach  starts  with  an  initial  design.  The  purpose  of  this  initial 

design  is  to  provide  preliminary  information  about  the  response  space.  Based  on  the 
initial  design,  LOO  errors  are  calculated  for  each  response  and  a  metamodel  built  for 
LOO  errors  corresponding  to  each  response.  As  mentioned  in  Chapter  3,  the  LOO 
error  calculation  used  in  this  dissertation  is  the  relative  LOO  error  defined  in  Chapter 
3.  A  constrained  multiobjective  optimization  problem  is  then  fonnulated  which  aims 
at  finding  points  in  the  design  space  that  simultaneously  have  the  maximum  value  of 
the  LOO  error  for  each  response  and  are  at  the  same  time  sufficiently  away  from 
existing  points  in  the  current  design,  D.  The  solution  to  this  multiobjective 
optimization  problem  is  a  set  of  Pareto  points  from  which  a  single  point  is  chosen  as 
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the  next  sample.  The  single  point  is  chosen  based  on  the  criterion  described  in  Section 
4.3.5.  The  newly  chosen  point  is  then  evaluated  for  its  true  response  and  is  added  to 
the  current  design  D.  The  procedure  is  repeated,  starting  with  the  LOO  error 
calculation  followed  by  choosing  the  next  sample,  until  a  stopping  criteria  is  reached 
(e.g.,  such  as  the  maximum  number  of  function  calls  is  reached).  The  details  of 
individual  steps  are  described  in  the  following  sections. 

4.3.2  Choice  of  Initial  Design 

As  mentioned  in  Chapter  3,  the  initial  design  is  crucial  to  the  success  of  an 
adaptive  DOE  strategy  and  the  resulting  metamodel.  In  the  proposed  approach,  it  is 
recommended  to  use  the  MED  method  to  generate  an  initial  design.  Again,  the  reason 
for  using  the  MED  method  is  to  have  a  good  representation  of  the  entire  input  space. 
The  number  of  points  in  the  initial  design  is  influenced  by  the  size  of  the  problem  as 
well  as  the  response  behavior.  Since  the  response  behavior  is  not  known  a  priori,  we 
use  ( 1  ()xd)  as  the  number  of  points  in  the  initial  design,  where  d  is  the  problem 
dimension  (i.e.,  number  of  input  variables). 

4.3.3  LOO  Error  Calculation 

In  Chapter  3  it  was  proposed  to  use  a  Kriging  based  LOO  error  prediction  model 
for  the  response.  In  case  of  multiple  responses,  we  extend  our  earlier  single-response 
approach  to  include  a  LOO  error  model  for  multiple  responses.  The  metamodeling 
technique  for  this  LOO  error  model  is  Kriging.  We  can  use  multiple  single  response 
Kriging  models,  a  multiresponse  Kriging  model  or  the  Cokriging  method  (see  Section 
4.2).  We  use  multiple  single-response  Kriging  metamodels  as  opposed  to  a 
multiresponse  Kriging  model.  In  other  words,  the  LOO  errors  for  each  response  are 
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calculated  separately  and  one  metamodel  for  LOO  errors  corresponding  to  each 
response  is  developed.  The  difference  is  that  when  individual  metamodels  are  used, 
the  correlation  parameters  are  calculated  for  each  response  separately.  In 
multiresponse  Kriging,  there  is  only  one  set  of  correlation  parameters  for  all  the 
responses,  resulting  in  less  accurate  metamodels.  The  Gaussian  correlation  discussed 
in  Chapter  2  is  used  for  Kriging,  but  other  correlations  can  be  used  and  in  some  cases 
may  be  required  due  to  the  nature  of  the  responses. 

4.3.4  Space-filling  Criterion 

The  space  filling  criterion  is  similar  to  the  one  proposed  in  Chapter  3.  The  metric 
S  is  defined  as  follows: 

ds(xf)  =  minfll  x  -x  |L  )  ,  V  x  e  D  a  (i  ^  /') 

(4.i) 

S  =  0.5  x  ma x(ds(xi)),  Vx;.  e  D 


4.3.5  Choosing  Next  Sample 

As  mentioned  in  Chapter  3,  the  problem  of  choosing  the  next  sample  point  is  a 
tradeoff  between  allocating  of  the  experiment  in  a  CAMM  region  (with  higher  LOO 
error  values  than  other  regions)  or  allocating  the  sample  in  a  sparsely  sampled  region 
of  the  space  (indicated  by  the  space  filling  metric).  In  order  to  extend  this  to  multiple 
responses,  we  choose  the  LOO  errors  for  different  responses  as  the  tradeoff  values 
subject  to  the  space  filling  metric  constraint.  This  is  mathematically  represented  as 
follows: 
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Note  that  the  above  optimization  problem  is  an  MOOP  and  thus  will  generally 
yield  a  set  of  Pareto  optimal  solutions.  The  challenge  then  is  to  select  one  or  more 
points  to  conduct  experiments.  There  are  different  methods  that  can  be  used  to  choose 
the  next  samples  as  follows: 

(a)  Choose  all  the  Pareto  points.  This  would  involve  sampling  a  large  number 
of  points  in  each  stage. 

(b)  Choose  points  from  the  Pareto  set  that  have  the  highest  value  of  the  LOO 
error  for  a  given  response,  implying  a  preference  or  affinity  towards  one 
particular  response. 

(c)  An  extension  of  (b)  is  to  choose  all  points  from  the  Pareto  set  that  have  the 
highest  predicted  LOO  error  values  for  each  of  the  responses,  resulting  in 
(at  most)  m  points  being  sampled  per  stage,  where  m  is  the  number  of 
responses. 

(d)  Use  a  more  elaborate  technique  such  as  level  of  dominance  approach 
described  in  Section  3.3.3. 

An  example  of  the  above  mentioned  tradeoff  between  the  LOO  errors  for  different 
responses  is  shown  in  Figure  4.1.  Figures  4.1a  and  4.1b  show  the  response  surfaces 
for  a  2-input  2-output  numerical  example.  The  initial  design  is  shown  in  Figure  4.1c. 
The  LOO  errors  are  calculated  for  all  points  in  the  design  space  and  using  a  non- 
dominated  sorting  algorithm,  a  Pareto  set  is  found  as  shown  in  Figure  4. Id.  As  seen 
from  Figure  4. Id,  clearly  there  is  a  tradeoff  between  the  LOO  errors  for  the  two 
responses.  The  corresponding  LOO  error  surfaces  for  the  two  responses  are  shown  in 
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Figures  4.1e  and  4. If.  We  observe  that  the  LOO  error  surfaces  are  nonlinear  and 
multimodal  in  nature. 

The  problem  in  Eq.  (4.2)  can  also  be  formulated  as  an  unconstrained  optimization 
problem  as  shown  in  Eq.  (4.3). 

V  =  argmax[^00,yj]  j  =  (4.3) 

X 

The  above  formulation  essentially  tries  to  maximize  2m  objectives,  where  m  is  the 
number  of  responses.  The  formulation  presented  in  Eq.  (4.2)  is  helpful  when  Kriging 
is  used  as  the  metamodeling  technique  and  when  individual  metamodels  are  built  for 
each  response,  because  the  Kriging  predictor  also  calculates  the  value  of  s  along  with 
the  predicted  response. 
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(e)  (f) 

Figure  4.1  LOO  error  tradeoff,  (a)  Response  surface  forji,  (b)  response  surface  for  y2, 
(c)  initial  design,  (d)ji  vs.j2  LOO  error  tradeoff,  (e)  LOO  error  surface  forji,  and  (f) 

LOO  error  surface  for 
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4.3.6  Design  Update 

The  new  point  or  points  selected  in  Section  4.3.5  are  evaluated  for  their  true 
responses  and  added  to  the  current  design  D. 

4.3.7  Stopping  Criteria 

The  stopping  criteria  in  this  case  can  be  one  or  more  of  the  following: 

(a)  Maximum  number  of  function  calls 

(b)  Acceptable  error  measure  (i.e.,  RMSE,  MAE)  or  performance  criteria  such  as 
MAS  (described  in  Chapter  2),  for  one  or  more  responses  based  on  an  existing 
test  sample  of  randomly  generated  points. 

(c)  Acceptable  value  for  the  norm  of  an  error  measure,  e.g.,  the  maximum  value 
of  MAE  amongst  all  responses,  based  on  a  test  sample. 

Generally  in  practical  applications,  the  maximum  number  of  available  function 
calls  is  used  as  the  stopping  criterion  and  the  same  is  used  in  this  dissertation. 

4.3.8  Step-by-Step  Description  of  the  MSFCVT  Approach 

Let  Kmax  be  the  maximum  number  of  available  function  calls.  The  different  steps 

are  as  follows: 

Step- 1 :  Generate  an  initial  design  D  with  predefined  number  of  samples  using  the 
MED  method. 

Step-2:  Calculate  LOO  errors  for  all  responses  for  all  the  points  in  the  initial 
design  and  build  the  Kriging  based  metamodel  for  the  LOO  errors. 

Step-3:  Solve  the  optimization  problem  discussed  in  Eq.  (4.2)  or  (4.3)  to  obtain 
the  next  sample(s). 
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Step-1  Step-2  Step-3 


Figure  4.2  Step-by-Step  flow  chart  for  MSFCVT  approach 


Step-4:  Check  if  the  optimization  problem  in  Step-3  is  solved.  If  the  optimization 
fails,  or  if  there  are  no  feasible  solutions  go  to  Step-7,  else  go  to  Step-5. 


Step-5:  Evaluate  the  true  response  for  the  new  point(s)  and  add  it  to  D. 


Step-6:  If  the  number  of  points  in  D  is  equal  to  the  available  number  of  function 
calls,  then  Stop,  otherwise  continue  to  Step-2. 


Step-7:  Compute  the  next  sample  using  either  the  MED  or  the  MSE  method  and 
proceed  to  Step-5. 
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4.4  Application  Examples 


In  this  section,  we  apply  the  proposed  approach  to  several  numerical  and 
engineering  test  examples  from  the  literature  and  from  the  test  suite  developed  as  a 
part  of  this  dissertation  to  evaluate  its  applicability  and  scalability  to  a  wide  range  of 
problems. 

The  metrics  of  MAS,  RMSE,  MAE,  RRMSE,  RMAE  and  the  Fol  are  calculated 
individually  for  each  response. 

4.4.1  Numerical  Example-1 

The  first  numerical  example  is  taken  from  Li  et  al.  (2006)  as  a  challenging  test 
problem  for  multiresponse  approximation.  We  note  that  the  response^  in  Eq.  (4.4)  is 
basically  the  product  of  x\  and  a  highly  nonlinear  function  (Farhang-Mehr  and 
Azarm,  2005)  which  was  described  in  Chapter  3.  This  makes  the  problem 
challenging.  The  problem  is  a  two-input  two-response  test  problem  with  a  closed- 
form  expression  and  is  shown  in  Eq.  (4.4). 

Vj  =  xlx2  +  0.5 

>2  =  xl{(l-e-2'^)  +  6x2e~7Xi  sin(10x2)- 
O.2r000(l2425)2  +  60 min 0,|  x2  -0.14 1  -0.08)2 
[ln(x2  +0.2)  +  1.5sin2(85x2)]} 

The  corresponding  response  surfaces  for  Vi  and  V2  are  shown  in  Figure  4.3. 

The  total  number  of  designs  was  fixed  at  20  for  a  fair  comparison  with  the  DMM 
(Li  et  al.,  2006;  Li,  2007)  design  and  RMSE  reported  in  Li  et  al.  (2006).  The 
MSFCVT  design  was  generated  starting  with  an  initial  design  of  10  points  using  the 
MED  method  and  then  adaptively  sampling  10  additional  points,  one  point  at  a  time. 
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A  set  of  2000  randomly  generated  points  (i.e.,  test  sample)  was  used  for  metamodel 
verification.  The  corresponding  designs  and  response  surface  for  the  MED  and  the 
DMM  (Li  et  al.,  2006)  method  are  shown  in  Figures  4.4  and  4.5  respectively.  The 
numerical  experiments  are  repeated  10  times  to  reduce  statistical  bias  and  thus  for 
each  metamodel  performance  metric,  we  have  a  range  of  values.  As  in  Chapter  3,  a 
box  plot  (McGill  et  al.,  1978)  is  used  to  concisely  display  the  results. 


Figure  4.3  Responses  for  numerical  example  -1  from  Eq.  (4.4),  (a)  Response  foryi,  and 

(b)  Response  fory2- 

From  Figure  4.4a.,  we  observe  that  the  points  generated  using  the  MED  method 
try  to  fill  the  space  rather  uniformly  while  the  MSFCVT  samples  (Figure  4.4b)  are 
more  concentrated  in  the  lower  left  of  the  X1-X2  space.  This  can  be  explained  based  on 
the  irregular  behavior  of  the  second  response  shown  in  Figure  4.3b.  The  response 
surfaces  generated  using  the  metamodels  are  also  shown  in  Figure  4.4.  We  observe 
that  the  response  surface  for  vi  closely  resembles  its  true  response  surface  shown  in 
Figure  4.3a. 
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Figure  4.4  MED  vs.  MSFCVT  approach,  resulting  designs  and  the  response  surfaces 
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In  case  of  V2,  the  response  surface  generated  using  the  MSFCVT  based 
metamodel  more  closely  resembled  the  true  response  surface.  These  observations  are 
further  supported  by  the  metamodel  performance  metrics  for  the  corresponding 
responses.  We  also  note  that  the  response  y  \  is  a  simple  function  and  as  such  after  20 
samples,  the  metamodel  performance  is  more  than  adequate  for  response  y  1,  as  the 
MAS  values  are  100  for  metamodels  resulting  from  both  the  MED  and  the  MSFCVT 
methods. 

Figure  4.5  shows  the  published  design  generated  using  the  DMM  (Li  et  ah,  2006) 
method  and  that  generated  by  the  MSFCVT  method.  The  design  from  the  MSFCVT 
method  is  the  same  as  from  Figure  4.4.  The  DMM  (Li,  et  ah,  2006)  method  on  the 
other  hand  appears  to  sample  points  away  from  the  center  of  the  design  space,  close 
towards  the  corners.  We  also  observe  that  the  metamodel  for  response  y\  is  near 
perfect  from  both  the  DMM  and  MSFCVT  methods.  As  for  response  V2,  again  it  can 
be  seen  that  the  generated  response  surface  from  MSFCVT  more  closely  resembles 
the  true  response  surface  than  the  one  from  DMM  (Li  et  ah,  2006)  method. 

The  MAS,  RRMSE  and  RMAE  comparisons  presented  in  Figures  4.6,  4.7  and  4.8 
in  the  form  of  box  plots.  Figure  4.6  shows  the  box  plots  for  the  MAS  obtained  using 
the  different  methods  for  the  two  responses.  The  RMSE,  RRMSE  and  RMAE 
obtained  from  the  3  methods  for  the  two  responses  Vi  and  \h  are  shown  in  Figures  4.7 
and  4.8  respectively.  Figure  4.9  shows  the  relative  errors  for  the  randomly  selected 
test  points  (for  1  run)  plotted  as  a  scatter  in  the  design  space. 

From  Figure  4.6,  it  can  be  seen  that  all  the  three  methods  viz.  MED,  DMM  and 
the  proposed  MSFCVT  yield  a  near  perfect  (i.e.,  MAS  =  100)  metamodel  for  Vi.  But 
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the  MSFCVT  approach  also  provides  better  RMSE  and  MAE  values  for  the  resulting 
metamodels  as  compared  to  DMM  and  MED  method  as  seen  from  Figure  4.7.  For  V2, 
based  on  the  median  MAS  values,  we  see  that  the  MSFCVT  approach  performs  better 
than  the  MED  and  the  DMM  (Li  et  ah,  2006)  approach. 

Based  on  median  values  (indicated  by  the  vertical  line  within  each  box)  from 
Figure  4.8,  we  observe  that  the  proposed  approach  outperforms  the  DMM  approach  in 
terms  of  RMSE,  RRMSE  and  RMAE,  but  is  inferior  to  MED  method  when  RRMSE 
and  RMAE  are  considered.  From  Figure  4.5,  we  see  that  the  MSFCVT  method 
allocates  more  points  in  the  lower  left  region  of  the  design  space.  The  resulting 
metamodel  from  MSFCVT  has  a  better  performance  in  this  region  which  is  verified 
by  the  corresponding  distribution  of  errors  as  a  function  of  test  points  as  shown  in 
Figure  4.9. 
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Figure  4.5  DMM  vs.  MSFCVT,  resulting  designs  and  metamodels 
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Figure  4.6  MAS  Comparison,  (a)  Response^  and  (b)  Response^ 
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Figure  4.7  RRMSE,  RMAE  and  RMSE  Fol  values  for  from  different  methods 
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Figure  4.8  RMSE  for  second  response^ 
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Figure  4.9  Error  distribution  forj2>  (a)  DMM  method,  and  (b)  MSFCVT  method 


The  RMSE  for  yi  as  reported  in  Li  et  al.  (2006)  for  the  DMM  method  is  0.11 
while  the  RMSE  obtained  from  the  proposed  approach  is  0.045,  thus  supporting  the 
conclusion  that  the  proposed  approach  outperforms  the  DMM  method  for  this 
example. 


4.4.2  Numerical  Example  -  2 

This  second  numerical  example  is  taken  from  Romero  et  al.  (2006)  to  allow  for  a 
fair  comparison  with  the  results  published  in  the  literature.  This  test  example  is  a 
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combination  of  two  single-response  test  examples  to  represent  a  multiresponse 
simulation.  It  has  two  inputs  and  two  outputs.  The  two  test  problems  are  Osio’s 
Function  and  Sasena’s  function,  designated  as  P9  and  P8  respectively  and  their  details 
are  provided  in  Appendix-A. 

Similar  to  the  first  example  in  the  earlier  section,  the  proposed  approach  is  applied 
to  this  test  problem.  The  approach  starts  out  with  9  initial  designs  and  a  total  of  19 
designs.  The  number  of  design  points  was  chosen  to  be  19  to  allow  comparison  with 
results  published  in  Romero  et  al.  (2006).  A  set  of  randomly  generated  2000  test 
points  is  used  for  metamodel  verification.  Changing  the  number  of  verification  points 
to  other  values  (i.e.,  1000,  5000)  resulted  in  a  similar  metamodel  performance. 

The  true  response  surface  for  these  two  examples  is  shown  in  Figure  4.10a  and 
4.10b  respectively.  Figure  4.11  shows  the  design  obtained  using  the  MED  and  the 
MSFCVT  methods  and  the  corresponding  response  surfaces  for  the  two  responses. 
From  Figures  4.1  la  and  4.1  lb,  we  observe  that  the  designs  generated  using  the  MED 
and  the  MSFCVT  methods  appear  to  be  space-filling  in  nature,  but  only  the  MED 
points  are  space-fdling.  The  corresponding  metamodel  generated  response  surfaces 
are  shown  in  Figures  4.1  lc  and  4.1  Id  respectively.  The  surface  generated  for  vi  using 
the  MSFCVT  method  resembles  more  closely  to  the  true  response  surface.  In  case  of 
V2,  both  methods  do  not  represent  the  surface  very  well. 
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Figure  4.10  Response  surfaces  for  Numerical  Example-2,  (a)  response  ji,  and  (b) 

response  y>2 

One  reason  for  V2  not  being  represented  very  well,  is  the  initial  design.  If  more 
number  of  points  were  to  be  added  to  the  initial  design,  the  predicted  response 
surfaces  would  be  closer  to  the  true  surfaces. 

The  corresponding  metamodel  performance  metrics  are  shown  in  Figures  4.12, 
4.13  and  4.14.  Figures  4.12  and  4.13  show  the  actual  RMSE  and  MAE  for  the  two 
methods  (i.e.,  MED  and  MSFCVT)  for  responses  vi  and  V2  respectively.  Figure  4.14 
shows  the  MAS  for  the  two  responses.  Based  on  Figures  4.12,  4.13  and  4.14,  we 
observe  that  the  MSFCVT  approach  outperforms  the  MED  method  for  response  vi, 
but  is  marginally  inferior  with  respect  to  V2- 
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(e)  (f) 

Figure  4.11  Numerical  Example-2,  (a)  MED  design,  (b)  MSFCVT  design,  (c)  MED 
response  ji,  (d)  MSFCVT  response yi,  (e)  MED  response^,  (f)  MSFCVT  response^ 
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Figure  4.12  Comparison  of  RMSE  and  MAE  for  response 

SFCVT  MAE 

SFCVT  RMSE 

MED  MAE 

MED  RMSE 


Figure  4.13  Comparison  of  RMSE  and  MAE  for  response^ 
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Figure  4.14  MAS  Comparison  forj]  andj2  from  MED  and  MSFCVT  approaches 
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Comparing  the  above  results  with  those  published  in  Romero  et  al.  (2006),  we 
observe  that  the  MSFCVT  approach  provides  a  better  MAE  than  the  three  approaches 
they  proposed. 


4.4.3  Application  to  Numerical  Test  Problems 

In  this  section,  a  suite  of  8  test  problems,  designated  at  PM1  to  PM8,  was 

complied  based  on  test  problems  from  the  literature.  The  details  of  this  multiresponse 
test  problem  suite  are  given  in  Appendix  A.  In  general,  several  single -response  test 
problems  are  combined  to  represent  a  multiresponse  computer  simulation.  A 
summary  of  the  numerical  experiments  that  were  conducted  is  shown  in  the  table 
below. 


Table  4.1  Test  Matrix  for  Multiresponse  Adaptive  DOE  Method 


Test 

Number  of 

Number  of 

Initial 

New 

Random 

Problem 

Inputs 

Outputs 

Design 

Samples 

Test  Points 

PM  1  -  PM2 

4 

2 

40 

40 

4000 

PM3  -  PM4 

8 

2 

80 

80 

8000 

PM5  -  PM6 

2 

4 

20 

20 

2000 

PM7 

4 

4 

40 

40 

4000 

PM  8 

8 

4 

80 

80 

8000 

Each  of  the  test  problems  (i.e.,  PM1  to  PM8)  has  multiple  responses.  We  use  the 
notation  PM-1  to  indicate  the  first  response  for  test  problem  PM1,  PM  1-2  the  second 
response  for  problem  PM  1  and  so  on. 

The  first  set  of  results  is  the  comparison  of  MAS  obtained  using  the  two  different 
selection  strategies  for  the  MSFCVT  method  using  Eq.  (4.2).  The  first  selection 
method  is  based  on  the  level  of  dominance  approach  and  selects  one  point  per  stage. 
The  second  selection  method  chooses  the  extreme  Pareto  points  at  each  stage  (i.e. 
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method  (c)  described  in  Section  4.3.5).  The  comparison  is  shown  in  Figure  4.15  for 
problems  PM-1  to  PM-4  and  in  Figure  4.16  for  problems  PM5-PM8. 
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Figure  4.15  MAS  comparison  for  MSFCVT  method,  choosing  1  point  vs.  2  points  per 

stage,  for  problems  PM1  -  PM4 
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Figure  4.16  MAS  comparison  for  MSFCVT  method,  1  point  vs.  2  points  per  stage,  for 

problems  PM5  -  PM8 
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From  Figure  4.15,  for  test  problems  PM1  to  PM4,  we  observe  that  there  is 
marginal  difference  between  the  MAS  using  the  two  methods.  On  the  other  hand,  in 
the  4-output  problems  (i.e.,  PM5-PM8)  there  is  significant  difference  in  MAS  for 
responses  PM5-4,  PM6-1,  PM6-3,  PM6-4,  PM7-1  and  PM7-3.  Overall  the  2  points 
per  stage  method  works  better  than  1  point  for  problems  PM6  and  PM7.  This  is 
mainly  because,  at  each  stage,  4  new  points  are  added  (as  opposed  to  one  point  in  the 
1  point  per  stage  approach),  which  also  increases  the  accuracy  of  the  intermediate 
LOO  error  models  for  each  of  the  responses.  This  results  in  a  better  overall 
metamodel  even,  when  the  total  number  function  calls  is  the  same.  The  RMSE  and 
MAE  results  were  comparable  for  the  two  responses  and  hence  are  not  shown  here. 

Based  on  these  test  problems,  the  proposed  MSFCVT  method  was  compared  with 
the  MED  method  and  the  results  are  discussed  next.  The  comparison  with  the  MED 
method  is  based  on  using  1  point  per  stage  for  the  MSFCVT  method.  The  MAS  for 
the  resulting  metamodels  are  shown  in  Figure  4.17.  The  proposed  approach  provides 
strictly  better  MAS  for  both  responses  for  2  out  of  4  test  problems.  The  corresponding 
Fol  results  for  RMSE,  MAE,  RRMSE  and  RMAE  Fol  are  shown  in  Figure  4.18. 
Except  for  the  second  response  of  problem  PM3,  the  proposed  approach  performs 
comparably  and  in  some  cases  even  better  (e.g.,  PM1)  than  the  MED  method. 
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Figure  4.17  MAS  for  MED  vs.  MSFCVT  method  for  2  response  case 
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Figure  4.18  MED  vs.  MSFCVT,  Fol  values  for  2  response  cases 
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Figure  4.19  MAS,  MED  vs.  MSFCVT  method  for  4  response  cases 


The  remaining  4  test  problems  (PM5-PM8)  are  problems  with  4  responses.  The 
results  for  these  are  shown  in  Figures  4.19  and  4.20.  From  Figure  4.19,  in  all  the  test 
problems,  except  for  PM6,  the  proposed  MSFCVT  approach  performs  comparable  to 
or  better  than  the  MED  method.  Overall,  based  on  the  MAS  values,  the  MSFCVT 
approach  gives  better  metamodels  for  10  out  of  16  responses.  For  PM6,  from 
Appendix- A,  we  note  that  the  4th  response  is  the  Busby  et  al.(2007)  problem  (PM  10) 
for  which  we  concluded  from  Chapter  3  that  the  implemented  SFCVT  approach  does 
not  work  well,  because  of  weak  metamodels  for  the  LOO  errors  in  intermediate 
stages.  Consequently,  the  MED  method  performs  significantly  better  in  terms  of 
MAS  as  seen  from  the  bars  for  PM6-4  in  Figure  4.19. 


From  Figure  4.20,  we  observe  that  the  MSFCVT  approach  gives  better  RMSE, 
MAE,  RRMSE  and  RMAE  for  7  out  of  16  responses  and  it  yields  better  RRMSE  and 
RMAE  in  1 1  out  of  16  responses. 
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Figure  4.20  MED  vs.  MSFCVT,  Fol  for  4  response  cases 

It  is  important  to  point  out  that  the  test  problems  PM5-PM8  used  the  single 
response  test  problems  from  Chapter  3  to  mimic  a  multiresponse  simulation  and  the 
default  implementation  of  the  MSFCVT  method  was  used.  Thus,  the  performance  of 
the  MSFCVT  method  when  compared  with  MED  method  will  be  similar  to  that 
observed  for  the  single  response  cases. 

4.4.4  Comparison  Between  Formulations 

In  this  section  we  briefly  present  and  compare  the  results  obtained  using  the  two 

different  formulations  given  by  Eq.  (4.2)  and  (4.3).  Note  that  there  are  m  objectives  in 
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the  formulation  presented  by  Eq.  (4.2),  where  as  there  as  2m  objectives  in  the 
formulation  presented  in  Eq.  (4.3),  where  m  is  the  number  of  responses.  Figure  4.21 
shows  the  MAS  obtained  by  applying  the  two  formulations  to  test  problems  PM1- 
PM4.  For  both  formulations,  the  extreme  Pareto  points  were  chosen  as  the  next 
samples.  For  the  first  response,  the  m  objective  formulation  seems  to  yield  better 
MAS  than  the  second  formulation.  Except  for  test  problem  PM-1,  the  two 
formulations  yield  close  results  for  the  other  3  test  problems.  It  should  be  noted  that 
even  though  extreme  Pareto  points  corresponding  to  each  objective  are  selected,  the 
number  of  unique  samples  is  not  the  same  as  the  number  of  objectives.  In  other 
words,  when  we  select  extreme  points  for  n  conflicting  objectives,  we  should  have  n 
solutions.  But  when  some  of  them  are  not  conflicting  (which  we  don’t  know  before 
hand),  then  the  extreme  points  for  two  or  more  objectives  coincide.  Thus  the  number 
of  unique  samples  chosen  is  not  the  same  as  the  number  of  problem  objectives. 
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Figure  4.21  Comparison  of  MAS  for  test  problems  PM1-PM4,  for  MSFCVT 
formulation  with  m  objectives  vs.  2 in  objectives 
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The  results  for  the  test  problems  with  four  responses,  i.e.,  test  problems  PM5- 
PM8,  are  shown  in  Figure  4.22.  The  MAS  values  are  mixed  (low  or  high  for  different 
responses)  for  test  problems  PM5-PM7,  while  for  problem  PM8  the  results  for  all 
responses  are  comparable. 
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Figure  4.22  Comparison  of  MAS  for  test  problems  PM5-PM8  for  MSFCVT  formulation 

with  m  objectives  vs.  2m  objectives 


Application  to  additional  test  problems  is  necessary  to  draw  any  generalized 
conclusions  with  regards  to  the  two  formulations. 

4.5  Summary 

In  this  chapter  a  new  adaptive  design  of  experiments  approach  for  multiresponse 
computer  simulations  is  presented.  The  approach  extends  the  concept  of  the  tradeoff 
between  a  space-filling  criteria  and  the  LOO  error  proposed  in  Chapter  3  to  multiple 
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responses.  In  this  case,  the  constrained  formulation  is  used  to  choose  the  next  sample, 
but  now  there  are  multiple,  generally  competing,  objectives.  One  or  more  points  can 
be  chosen  as  next  samples.  The  new  sample  chosen  needs  to  a  strike  a  balance 
between  improving  the  metamodels  for  all  the  different  responses  simultaneously  vs. 
improving  metamodels  for  individual  responses.  Alternatively,  all  the  Pareto  points 
could  be  chosen  as  the  next  samples.  While  this  approach  will  utilize  more  function 
calls  per  stage,  it  will  give  a  better  representation  of  the  response  space  in  the  early 
stages.  The  drawback  of  choosing  all  Pareto  points  at  each  stage  is  the  possibility  that 
the  new  points  will  be  allocated  much  closer  to  existing  samples.  The  third  option  to 
select  new  samples  from  the  obtained  Pareto  set  is  simply  to  choose  the  extreme 
points.  This  will  result  in  the  selection  of  multiple  samples  per  stage  with  the  number 
of  samples  equal  to  the  number  of  responses. 

The  proposed  approach  is  compared  with  two  multiresponse  numerical  examples 
from  the  literature.  The  first  test  problem  is  from  Li  et  al.  (2006)  and  it  was  found  that 
the  proposed  approach  outperforms  the  DMM  method  proposed  by  the  Li  et  al.  when 
applied  to  their  particular  test  problem.  The  approach  was  also  compared  against 
published  results  for  three  other  multiresponse  adaptive  DOE  approaches  proposed  by 
and  applied  to  one  test  example  in  Romero  et  al.  (2006).  Only  the  MAE  data  was 
available  and  it  was  found  that  the  MSFCVT  approach  generates  design  which  result 
in  better  metamodels  than  those  published  in  Romero  et  al.  (2006).  The  proposed 
approach  was  compared  to  the  MED  method  for  a  suite  of  8  test  problems.  It  was 
found  that  the  proposed  MSFCVT  method  outperforms  the  MED  method  in  majority 
of  the  test  problem  responses. 
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In  the  next  chapter,  we  present  an  approximation  based  framework  for  the  multi¬ 
level  performance  evaluation  of  novel  heat  exchangers  and  apply  the  single  and  the 
multiresponse  adaptive  DOE  approaches  to  the  heat  exchanger  design  problem. 


131 


Chapter  5  Approximation  Based  Optimization  Framework  for 
Design  of  Novel  Heat  Exchangers 

5.1  Introduction 

In  this  chapter,  a  new  framework  for  the  multi-level  performance  evaluation  and 
design  optimization  of  air-cooled  heat  exchangers  is  presented.  The  framework  uses 
metamodels  for  CFD  simulations,  using  the  adaptive  DOE  approach  developed  in 
Chapter  3  for  design  optimization  of  novel  air-to-refrigerant  heat  exchangers. 

Heat  Exchangers  (HX)  are  used  to  transfer  energy  from  one  working  fluid  or 
medium  to  another  using  a  temperature  gradient  as  the  driver.  The  working  fluids  can 
be  air,  water,  glycols  or  other  refrigerants  such  as  R134a.  Air-to-refrigerant  heat 
exchangers  are  widely  used  today  in  refrigeration  and  Air-Conditioning  (AC) 
applications.  The  use  of  these  heat  exchangers  include:  (a)  Indoor  and  outdoor  coils 
for  residential  and  commercial  buildings  such  as  houses,  offices,  super  markets  (b) 
the  condenser,  radiator  and  the  evaporator  in  an  automobile  air  cooling  system  (c)  the 
coil  on  the  back  or  bottom  of  a  household  refrigerator,  (d)  charge  air  coolers  and 
intercoolers  in  vehicles  etc.  Most  of  the  cooling  equipment  is  based  on  the 
thennodynamic  cycle  termed  as  the  vapor  compression  cycle  or  some  variation  of  it. 
In  general,  any  vapor  compression  cycle  based  equipment  requires  a  minimum  of  two 
heat  exchangers.  In  a  basic  vapor  compression  cycle  there  are  four  components  viz. 
(a)  compressor,  (b)  high-pressure  heat  exchanger,  termed  as  condenser  (c)  an 
expansion  device  and  (d)  a  low  pressure  heat  exchanger  termed  as  evaporator. 
Together,  the  two  heat  exchangers  in  an  AC  unit  account  for  the  more  than  50%  for 
the  initial  cost  (manufacturing  cost,  also  known  as  the  first  cost)  and  significantly 
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affect  the  performance  and  hence  the  energy  consumed  by  the  unit.  In  mobile 
applications,  the  weight  of  these  heat  exchangers  also  plays  an  important  role  in 
addition  to  the  first  cost,  since  it  affects  the  fuel  efficiency  of  the  vehicle. 

In  the  light  of  energy  efficiency  and  energy  conservation  concerns,  there  is  a 
growing  need  for  more  compact,  high  performance  heat  exchangers  and  ever  present 
requirement  to  minimize  cost.  In  the  last  decade  a  significant  increase  in  compactness 
has  been  achieved  in  air-to-refrigerant  heat  exchangers  by  replacing  the  traditional 
round-tube  plate  fin  heat  exchangers  with  microchannel  heat  exchangers. 

The  optimization  framework  proposed  in  this  chapter  aims  at  developing  a  novel 
tool  for  multi-level  design  optimization  of  new  heat  exchangers,  termed  as  the  next 
generation  heat  exchangers  (NGHX).  The  novel  aspects  of  this  approach  are  (a)  to  use 
simplified  2D  CFD  to  evaluate  performance  of  novel  heat  exchanger  elements  and  (b) 
to  develop  and  use  an  s-NTU  based  heat  exchanger  model  to  evaluate  component 
level  performance  of  the  novel  heat  exchanger  elements  and  (c)  to  use  the  adaptive 
DOE  approach  developed  in  Chapter  3  to  optimize  the  design  of  these  novel  heat 
exchangers.  Together,  these  novel  aspects  will  lead  to  better  heat  exchanger  designs 
which  are  more  compact  than  the  current  state  of  the  art,  as  will  be  shown  in  this 
chapter.  More  explanation  is  needed  to  highlight  the  multi-level  performance 
evaluation  aspect  of  this  research.  In  (a)  above,  the  performance  of  a  small  heat 
exchanger  element  is  evaluated  using  2D  CFD.  These  elements  are  then  “assembled” 
together  to  form  real  life  heat  exchangers  and  a  s-NTU  based  model  is  used  to 
evaluate  the  component  level  (i.e.,  heat  exchanger  as  a  whole)  performance  of  these 
heat  exchangers. 
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The  remainder  of  this  chapter  is  organized  as  follows:  Section  5.2  provides  some 
background  information  on  heat  exchangers  and  their  optimization.  Section  5.3  gives 
a  summary  of  literature  which  focuses  on  the  use  of  approximation  assisted 
optimization  for  heat  exchanger  design.  Section  5.4  provides  details  about  a  new  heat 
exchanger  geometry  that  will  be  analyzed  in  this  chapter.  Section  5.5  summarizes  a 
novel  approach  for  the  multi-level  performance  evaluation  of  the  new  heat  exchanger. 
Section  5.6  applies  the  adaptive  DOE  methods  developed  in  Chapters  3  to  the  heat 
exchanger  design  problem.  Section  5.7  provides  results  from  the  approximation 
assisted  optimization  for  the  NGHX.  Finally  Section  5.8  summarizes  this  chapter. 

5.2  Background  information 

In  this  section,  background  information  on  conventional  air-to-refrigerant  heat 
exchangers  in  presented.  An  introduction  to  an  s-NTU  based  simulation  tool  (i.e., 
CoilDesigner)  for  performance  evaluation  of  tube-fin  and  microchannel  heat 
exchangers  is  also  presented. 

5.2.1  Tube-Fin  Heat  Exchangers 

Tube-fin  heat  exchangers  (ARI,  2001)  generally  refer  to  round  tube  bundles  with 
fins.  The  fins  help  improve  the  air-side  heat  exchange  area.  A  schematic  of  a  tube-fin 
heat  exchanger  is  shown  in  Figure  5.1. 
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Figure  5.1  Schematic  of  a  tube-fin  heat  exchanger 


In  tube-fin  heat  exchangers,  tubes  are  arranged  in  a  bundle  across  fins.  The  tubes 
are  arranged  from  top  to  bottom  in  vertical  banks,  normal  to  the  direction  of  the  air 
flow.  The  refrigerant  (e.g.,  water)  flows  through  the  tubes  and  heats  or  cools  the  air 
based  on  the  application.  The  fins  mounted  on  the  tubes  are  used  to  extend  the  heat 
transfer  area  since  the  air  side  heat  transfer  coefficient  is  much  less  than  the 
refrigerant  side  heat  transfer  coefficient.  For  simplicity,  the  details  of  fin  patterns  are 
not  described  here. 

A  general  observation  from  Figure  5.1  is  that  there  is  a  significant  obstruction  to 
the  air  flow  due  to  the  tubes.  The  nominal  size  of  this  obstruction  in  the  air  path  is  the 
same  as  the  tube  diameter.  Any  improvement  in  this  air  flow  path  should  increase  the 
performance  of  the  heat  exchanger.  Using  this  concept  as  one  of  the  foundations,  the 
new  heat  exchanger  technology  conceived  after  tube-fin  heat  exchangers  were  heat 
exchangers  with  flat  tubes.  Flat  tubes  later  enabled  the  microchannel  heat  exchanger 
which  is  described  next. 
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5.2.2  MicroChannel  Heat  Exchangers 

A  schematic  of  a  microchannel  heat  exchanger  is  shown  in  Figure  5.2.  A 
microchannel  heat  exchanger  (Kays  and  London,  1984;  Kandlikar  et  ah,  2006) 
comprises  of  corrugated  fin  sheets  sandwiched  between  flat  tubes.  Each  tube  has  tiny 
(~  1  mm)  circular  or  rectangular  ports  through  which  the  refrigerant  flows. 
MicroChannel  heat  exchangers  are  generally  fed  using  an  inlet  and  an  outlet  header. 
The  division  of  each  tube  into  multiple  parallel  flow  channels  greatly  reduces  the 
pressure  drop  on  the  refrigerant  side.  Several  novel  fin  designs  (e.g.,  louver  fins)  are 
available  which  drastically  improve  the  air-side  performance  as  well.  The  first  view 
in  Figure  5.2  is  the  coil  face  area,  which  is  the  area  of  the  coil  which  ‘sees’  the  air 
first. 


Figure  5.2  Schematic  of  a  microchannel  heat  exchanger 
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5.2.3  CoilDesigner 

CoilDesigner  (Jiang  et  al.,  2006)  is  a  control  volume  based  simulation  tool  that 
can  simulate  the  performance  of  air-to-refrigerant  and  refrigerant-to-refrigerant  heat 
exchangers.  The  underlying  solver  allows  for  the  decomposition  of  individual  tubes  in 
a  heat  exchanger  into  smaller  heat  exchanger  elements  termed  as  ‘segments’. 
CoilDesigner  internally  uses  this  decomposition  along  with  the  e-NTU  method  (Shah 
and  Sekulic,  2003)  of  heat  transfer  calculations  which  helps  to  account  for  the  drastic 
changes  in  transport  properties  (density  and  viscosity)  during  evaporation  and 
condensation  processes  resulting  in  accurate  prediction  of  heat  exchanger 
performance.  The  methodology  behind  CoilDesigner  is  described  in  more  details 
along  with  experimental  validation  in  Jiang  et  al.  (2006),  Schwentker  et  al.  (2005)  and 
more  recently  in  Singh  et  al.  (2008).  Optimization  of  tube-fin  coils  using  MOGA  with 
CoilDesigner  has  been  described  and  demonstrated  in  Aute  et  al.  (2004).  It  is 
important  to  note  that  CoilDesigner  uses  correlations  for  air  and  refrigerant  side  heat 
transfer  and  pressure  drop  calculations.  This  makes  CoilDesigner  very  flexible, 
because  once  a  correlation  is  available  for  a  given  tube/fin  geometry,  CoilDesigner 
can  be  used  to  simulate  the  coil  performance  using  the  particular  tube/fin  geometry. 
Correlations  are  generally  based  on  experimental  data  sets,  but  in  cases  where 
experimental  data  is  not  available  (e.g.,  for  new  heat  exchangers),  the  heat  transfer 
and  pressure  drop  characteristics  can  be  obtained  using  a  detailed  CFD  simulation. 

In  order  to  evaluate  performance  of  heat  exchangers,  CoilDesigner  requires 
detailed  geometrical  and  design  information  such  as  tube  diameters,  thickness,  fin 
thickness,  tube  horizontal  and  vertical  spacing,  tube  length,  number  of  parallel  tubes, 
etc.  The  results  predicted  by  CoilDesigner  include  the  overall  heat  load,  refrigerant 
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side  heat  transfer  coefficients  (HTC),  refrigerant  side  pressure  drop  (RDP),  outlet 
refrigerant  and  air  conditions,  the  volume  of  heat  exchanger,  material  weight  etc. 

5.3  Literature  Review 

Conventional  heat  exchangers  such  as  tube-fin  and  microchannels  have  been  very 
well  studied  and  as  such  there  is  plenty  of  literature  which  focuses  on  the 
optimization  of  these  heat  exchangers.  In  terms  of  design  optimization  of  novel  heat 
exchangers,  the  literature  focuses  only  on  the  optimization  of  a  small  element  or  a 
segment  of  the  heat  exchanger.  The  published  research  can  also  be  classified  to  fall  in 
the  category  of  shape  optimization,  because  the  research  mainly  focuses  on  the 
optimization  of  a  fin  or  a  tube  shape.  Heat  Exchanger  optimization  can  be  carried  out 
at  the  element  level  (i.e.,  small  segment  of  heat  exchanger  is  analyzed),  at  the 
component  level  (i.e.,  an  entire  heat  exchanger  is  optimized  for  a  given  performance) 
and  at  the  vapor  compression  system  level  (i.e.,  heat  exchanger  geometries  are 
optimized  as  a  part  of  system  optimization).  While  there  is  ample  literature  available 
on  optimization  of  heat  exchangers  based  on  thennal-hydraulic  models,  the  literature 
review  presented  here  is  limited  to  the  application  of  approximation  assisted 
optimization  with  evolutionary  algorithms  to  the  design  optimization  of  novel  heat 
exchanger  geometries. 

The  use  of  CFD  based  optimization  is  fairly  common  in  the  aerospace  field.  Early 
examples  of  CFD  based  optimization  in  the  aerospace  field  include  the  classical 
aerofoil  optimization  (Falco,  1997;  Makinen,  et  ah,  1998)  problem.  CFD  based 
optimization  is  yet  to  spread  in  the  heating,  ventilation  and  the  air-conditioning  field, 
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and  remains  an  area  of  active  research.  Hilbert  et  al.  (2006)  use  MOGA  and  CFD  to 
optimize  the  shape  of  a  tube-fin  heat  exchanger  tube  for  maximum  heat  exchange 
while  obtaining  minimum  pressure  loss.  The  candidate  shapes  were  represented  via 
parameterized  NURBS  (Farin,  1999,  2002)  and  their  approach  used  CFD  coupled 
with  optimization  with  no  approximation  involved.  The  resulting  tube  designs,  similar 
to  aerofoil  shape,  were  highly  unconventional  yet  intuitive. 

Langer  et  al.  (2002)  propose  an  approach  for  shape  optimization  using 
parameterized  CAD  geometry  and  evolutionary  algorithms.  They  use  MATLAB  as 
the  integration  platform.  The  optimization  algorithm  generates  new  design  values  for 
the  inputs  required  for  a  parameterized  geometry  which  is  then  built  by  a  CAD 
program  followed  by  FEM  analysis  for  performance  evaluation.  While  the  integration 
approach  is  novel,  there  was  no  approximation  used  and  actual  FEM  runs  were 
carried  out  using  a  parallel  cluster  of  computers. 

Perez-Segarra  et  al.  (2002)  conducted  CFD  studies  on  conventional  tube-fin  heat 
exchangers  to  assess  the  reliability  of  the  numerical  solutions  obtained  from  the  CFD 
codes. 

Kanaris  et  al.  (2004)  propose  the  use  CFD  for  the  design  of  novel  flow  channels 
for  plate  heat  exchangers.  In  their  study,  the  CFD  results  are  first  validated  against 
experimental  data  for  known  channel  geometries.  Then,  new  geometries  are  explored 
in  a  parametric  fashion  in  order  to  find  better  flow  channels.  A  more  recent 
application  is  presented  in  Park  et  al.  (2006),  in  which  Kriging  and  CFD  is  used  for 
optimization  of  a  heat  sink. 
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Based  on  the  available  literature,  it  can  be  concluded  that  detailed  CFD  analysis 
and  optimization  is  carried  out  only  on  a  case-by-case  basis,  i.e.,  for  fixed  geometries. 
Moreover,  the  CFD  analysis  was  carried  out  at  the  element  level  and  no  studies  have 
been  conducted  to  evaluate  the  scalability  of  the  performance  at  the  component  (i.e., 
heat  exchanger)  level.  There  is  no  systematic  method  or  framework  for  the  analysis  of 
novel  heat  exchangers.  Morever,  approximation  assisted  optimization  is  seldom  used 
for  CFD  based  heat  exchanger  optimization. 

5.4  Novel  Heat  Exchanger  Concepts 

Heat  exchanger  design  optimization  is  in  general  a  very  broad  problem  and  it  is 
important  to  narrow  down  the  problem  definition.  As  a  starting  point  for  the  search  of 
NGHX,  it  is  important  to  develop  new  heat  exchanger  concepts.  These  concepts  can 
then  be  analyzed  further  using  CFD  and  approximation  assisted  optimization. 

After  a  comprehensive  literature  review,  a  workshop  to  brainstorm  for  new  heat 
exchanger  configurations  was  conducted  (Radermacher  et  ah,  2007).  The  result  of 
this  workshop  were  a  of  set  25  novel  heat  exchanger  configurations,  each  with  its 
own  geometry  variations.  Some  of  the  sample  geometries  are  shown  in  Figure  5.3 
below.  The  arrows  in  Figures  5.3a  and  5.3b  indicate  direction  of  the  air  flow. 

5.4.1  Potential  Geometries 

In  order  to  finalize  potential  geometries  for  further  analysis,  several  designs  were 
brainstormed  based  on  creativity.  Without  going  into  the  details  of  these  designs, 
some  basic  designs  were  chosen  for  detailed  analysis.  These  detailed  designs  are 
derived  from  the  concept  of  microchannel  heat  exchangers  and  are  based  on  the 
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notion  of  augmenting  the  heat  transfer  area  per  unit  external  volume  of  heat 
exchangers. 
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Figure  5.3  Concept  heat  exchanger  geometries,  Courtesy  :  NGHX  Project 
(Radermacher  et  al.,  2007),  (a)  and  (b)  multiport  staggered  configurations,  (c)  circular 

ports  with  fin  sheet 


5.4.2  Chosen  Geometry 

The  heat  exchanger  geometry  chosen  for  further  analysis  is  shown  in  Figure  5.4. 
This  geometry  has  several  round  tubes  connected  by  thin  sheets  that  act  as  a  tin.  Note 
that  compared  to  traditional  heat  exchangers  such  as  tube-fin  (reviewed  in  Section 
5.2.1)  and  microchannels  (reviewed  in  Section  5.2.2),  this  geometry  does  not  have 
any  fins,  though  the  metal  webs  connecting  the  individual  tubes  act  as  fins. 
Preliminary  analysis  (Radermacher  and  Azarm,  2007)  showed  that  this  geometry  is 
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capable  of  packing  more  heat  transfer  area  per  unit  external  volume  of  the  heat 
exchanger  per  unit  material  weight  of  the  heat  exchanger  as  compared  to  traditional 
microchannel  heat  exchanger,  when  port  diameter  is  less  than  0.4mm.  The  defining 
dimensions  of  the  geometry  are  inner  diameter  (Din),  horizontal  spacing  (Hs),  vertical 
spacing  (Vs)  and  the  length  (L)  and  depth  (w). 


(a) 


Figure  5.4  Proposed  NGHX  Geometry,  (a)  2-D  view,  (b)  3-D  view 
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Note  that  for  simplicity  we  choose  the  geometry  with  circular  flow  channels  (i.e., 
ports).  The  same  geometry  can  be  extended  to  include  any  other  arbitrary  port  shapes. 

5.4.3  Performance  Evaluation  Criteria 

Several  criteria  for  the  comparison  of  heat  exchangers  have  been  discussed  in  the 
literature  (Shah  and  Sekulic,  2003;  Webb  and  Kim,  2005).  These  criteria  are  mainly 
based  on  thermodynamic  characteristics  of  the  heat  exchanger.  The  particular 
criterion  of  interest  is  the  VG  criteria  (Webb  and  Kim,  2005)  in  which  heat 
exchangers  are  sized  for  a  required  thermal  heat  load  with  a  specified  flow  rate.  In 
this  case  study,  as  the  title  implies,  we  want  to  find  the  next  generation  of  air-cooled 
heat  exchangers  and  as  such  a  more  practical  comparison  criteria  was  chosen.  The 
two  criteria  that  were  chosen  were  the  total  heat  exchanger  volume  and  the  air-side 
pressure  drop  for  the  heat  exchanger  for  a  given  heat  duty.  A  1  kW  air-to-water 
microchannel  heat  exchanger  was  used  as  a  baseline  heat  exchanger  for  comparison. 
The  baseline  microchannel  was  designed  with  CoilDesigner  using  the  best  performing 
tube  and  fin  surface  configurations.  The  inlet  air  and  water  conditions  (pressure, 
temperature)  and  the  flow  rate  are  fixed. 

5.5  Performance  Evaluation  and  Design  Optimization  Approach 

In  this  section,  a  new  approach  for  multi-level  performance  evaluation  of  novel 
air-cooled  heat  exchangers  is  presented.  First,  the  underlying  assumptions  are  listed, 
followed  by  the  details  of  the  proposed  approach. 
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5.5.1  Assumptions 

The  assumptions  underlying  this  analysis  are  as  follows: 

(a)  The  heat  exchangers  can  be  represented  as  cross-flow  air-to-refrigerant 
heat  exchangers  as  shown  in  Figure  5.4.  A  single  tube  can  have  multiple 
ports  in  the  air-flow  direction  and  a  heat  exchanger  is  comprised  of 
multiple  tubes. 

(b)  The  working  fluids  are  air,  on  the  outside,  and  water  in  the  tubes.  While 
the  proposed  approach  is  applicable  for  other  fluids  such  as  R134a  instead 
of  water,  the  corresponding  analysis  is  much  more  complex  and  will 
include  two-phase  flow  for  which  heat  transfer  and  pressure  drop 
correlations  are  still  in  development.  Furthermore,  the  CFD  simulation  of 
two-phase  flows  is  still  an  open  research  problem.  Hence  water  is  used  to 
demonstrate  the  concept. 

(c)  The  performance  analysis  is  carried  out  for  heat  exchangers  comprising  of 
single  banks  only,  i.e.  one  set  of  tubes  from  top  to  bottom. 

(d)  The  performance  evaluation  criteria  are  the  ones  described  in  Section 
5.4.3.  Furthermore,  the  new  designs  are  compared  to  a  baseline 
microchannel  coil  that  can  deliver  the  same  performance. 

(e)  The  header  volume  and  header  performance  (i.e.,  flow  mal-distribution)  is 
ignored  in  this  analysis  and  as  such  while  comparing  the  coil  face  area, 
coil  volume  and  material  volume,  the  headers  are  ignored. 
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(f)  A  CFD  package  (Fluent,  2005)  will  be  used  for  any  performance 
calculations  that  require  CFD  analysis. 

(g)  In  the  current  analysis,  the  tube’s  inner  wall  is  assumed  to  have  a  constant 
temperature  for  the  length  of  the  segment  that  is  being  analyzed  using 
CFD.  In  an  actual  heat  exchanger,  each  segment  of  the  heat  exchanger 
will  have  a  different  tube-wall  temperature.  Numerical  experiments  were 
conducted  and  verified  that  this  assumption  does  not  have  a  significant 
effect  on  the  calculated  air-side  heat  transfer  performance  for  the  given 
range  of  input  conditions. 

5.5.2  Hybrid  Performance  Evaluation  Approach 

In  order  to  evaluate  the  performance  of  a  novel  HX  design,  there  are  two  main 

methods:  (a)  build  a  prototype  and  test  it  in  the  laboratory  and  (b)  model  the  HX  with 
a  detailed  simulation  tool  such  as  finite  element  method  (FEM)  or  computational  fluid 
dynamics  (CFD).  CFD  in  turn  can  be  carried  out  as  a  two  dimensional  (2D)  or  three 
dimensional  (3D)  analysis  based  on  the  geometry.  The  method  of  building  a 
prototype  for  each  design  is  not  feasible  due  to  resource  and  time  constraints,  which 
leaves  us  the  CFD  approach. 

The  evaluation  of  a  novel  heat  exchanger  with  the  current  simulation  technology, 
i.e.,  with  3D  CFD  requires  significant  computational  effort  (e.g.,  several  days  to 
weeks)  to  evaluate  a  single  assembled  heat  exchanger.  3D  Analysis  will  also  add  an 
additional  dimension  i.e.,  tube  length,  to  the  problem. 

The  proposed  hybrid  approach  aims  at  reducing  the  time  required  to  analyze  the 
performance  of  novel  heat  exchanger  elements  by  an  order  of  magnitude.  This 
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approach  uses  2D  CFD  for  a  small  number  of  cases  to  evaluate  the  heat  transfer 
coefficients  and  pressure  drops. 

The  high  level  steps  in  the  approach  can  be  summarized  as  follows: 

(a)  Evaluate  the  air-side  perfonnance  of  a  small  periodic  (repeating)  heat 
exchanger  element  using  2D  analysis.  HX  performance  requires  simultaneous 
analysis  of  fluid  flow  inside  the  tubes  and  the  air- flow  over  the  tubes.  Here  the 
assumption  is  that  the  tube-side  performance  is  readily  available,  as  is  the  case  for 
most  applications.  If  the  tube-side  performance  is  not  readily  available,  a  3D  analysis 
on  the  tube-side  can  be  conducted  to  get  the  performance.  Single  phase  fluid  flow  in 
microchannels  has  been  widely  investigated  (Morini,  2006;  Abdelaziz  et  al.,  2008) 
and  it  was  verified  that  the  conventional  correlations  are  still  valid  for  heat  transfer 
and  pressure  drop  calculations. 

(b)  Use  the  performance  evaluated  in  step  (a)  in  an  e-NTU  based  tool  such  as 
CoilDesigner  to  assemble  heat  exchangers  and  evaluate  its  performance.  The  new 
heat  exchanger  element  parameters  can  then  be  optimized  based  on  the  performance 
of  the  assembled  heat  exchanger. 

Thus,  we  evaluate  the  performance  of  the  novel  HX  element  at  the  finite  level  and 
use  the  finite  level  performance  information  along  with  CoilDesigner  to  model  the 
performance  of  an  entire  HX,  i.e.,  macro  level.  This  points  out  the  multi-level  aspect 
of  the  proposed  approach. 

The  air-side  performance,  i.e.,  air  heat  transfer  coefficient  (HTC)  and  air  pressure 
drop  (DP),  is  evaluated  for  a  small  segment  of  a  heat  exchanger  for  a  fixed  tube  inner- 
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wall  temperature.  It  has  been  shown  that  the  fixed  tube  inner-wall  temperature 
assumption  does  not  have  a  significant  impact  on  the  calculated  air-side  heat  transfer 
coefficient  values.  The  use  of  fixed  tube  inner  wall  temperature  also  implies  that  we 
do  not  need  to  account  for  the  flow  inside  the  tubes  and  this  reduces  the 
computational  effort  greatly. 

A  CoilDesigner  model  is  developed  that  can  evaluate  the  performance  of  the  new 
heat  exchangers.  This  new  model  can  use  the  air-side  performance  evaluated  earlier 
and  calculate  the  overall  performance  of  a  heat  exchanger.  In  a  separate  study 
(Abdelaziz  et  ah,  2008)  ,  it  was  shown  that  the  conventional  correlations  for  water 
side  heat  transfer  and  pressure  drop  calculations  are  valid  for  the  novel  heat 
exchanger  geometry  that  was  analyzed  in  this  chapter. 

5.5.3  Geometry  Representation  and  Parameterization 

Typically  all  CFD  tools  follow  a  three  step  procedure  (Versteeg  and 

Malalasekera,  1995)  as  follows:  (a)  Preprocessing  (geometry  development,  meshing 
and  boundary  conditions  setup),  (b)  Processing  (solving  the  governing  equations)  and 
(c)  Post  processing  (result  analysis  and  plotting).  In  the  geometry  development  phase, 
a  2D  model  (for  2D  analysis)  or  a  solid  model  (for  3D  analysis)  is  developed.  For  fast 
performance  evaluation  in  batch  mode,  the  geometry  should  be  parameterized  such 
that  by  changing  a  set  of  key  parameters,  the  shape  of  the  geometry  can  be  tweaked. 
An  example  is  shown  in  Figure  5.5.  By  changing  the  dimensions  a  and  b,  the  shape 
can  be  changed  from  a  circle  to  an  ellipse  and  vice-versa. 
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(a)  (b) 

Figure  5.5  Simple  shape  representation  and  manipulation 

The  analyses  of  these  geometries  require  perimeter  and  cross  section  area 
calculations.  While  such  calculations  are  trivial  for  simple  shapes  like  a  circle  or  a 
square  (as  shown  in  Figure  5.6a),  these  calculations  can  be  very  challenging  and  time 
consuming  for  more  complex  shapes  (the  ones  shown  in  Figure  5.6b).  Therefore  a 
sophisticated  method  of  representation  is  required  for  developing  geometries. 

It  is  proposed  to  use  non-unifonn  rational  B-splines  (NURBS)  to  represent  the 
different  geometries.  More  details  about  NURBS  can  be  found  in  Farin  (1999,  2002). 
The  advantages  of  NURBS  are  as  follows,  (a)  they  can  represent  practically  any  type 
of  curve,  i.e.,  smooth  or  with  kinks  (sharp  edges),  (b)  they  are  easy  to  manipulate  i.e. 
the  parameters  in  the  NURBS  equations  can  be  changed  to  change  the  shape  of  the 
curve  they  represent,  and  (c)  they  can  be  used  in  CFD  geometry  development  easily, 
as  many  popular  CFD  packages  support  NURBS. 

The  use  of  NURBS  will  allow  us  to  represent  the  geometry  in  a  parametric  form 
and  could  potentially  allow  the  design  optimization  at  the  port  level.  In  other  words, 
the  optimal  port  shape  could  be  an  aero-foil  or  an  elongated  triangle. 
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Even  though  it  is  proposed  to  use  NURBS  for  representing  shapes,  one  can  still 
use  simple  shapes  such  as  circles,  squares  etc.  if  the  HX  in  consideration  warrants  it. 
In  the  heat  exchanger  analysis  presented  in  this  chapter,  the  heat  exchanger  has 
circular  ports  and  they  are  specified  using  a  center  point  along  with  inner  and  outer 
radii. 

This  geometry  aspect  is  worth  mentioning,  since  when  it  is  coupled  with  the 
approach  proposed  in  this  chapter,  it  opens  up  a  new  avenue  for  rapid  and  streamlined 
evaluation  of  novel  heat  exchanger  geometries. 


(a) 


O 

o 


(b) 

Figure  5.6  Different  port  shapes,  (a)  simple  ports  and  (b)  complex  ports 
5.5.4  CFD  Analysis 

For  rapid  CFD  evaluation  of  novel  geometries,  the  process  of  geometry 
generation,  meshing  and  simulation  need  to  be  automated.  An  automated  tool 
(Radermacher  et  ah,  2007),  termed  as  Parameterized  Parallel  Computational  Fluid 
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Dynamics  (PPCFD)  was  used  to  carry  out  CFD  analysis  automatically  in  batch  mode. 
The  flow  chart  for  the  operation  of  this  PPCFD  tool  is  shown  in  Figure  5.7. 


Figure  5.7  Automated  CFD  evaluations 


Conventionally,  the  geometry  generation  and  the  mesh  generation  steps  are 
performed  manually  and  are  the  most  time  consuming.  In  the  new  tool  these  steps 
have  been  automated  by  a  computer  program  that  reads  input  geometry  parameters 
and  generates  these  files.  These  files  then  serve  as  input  to  the  CFD  solver.  The 
PPCFD  tool  itself  does  not  require  any  user  intervention  and  can  generate  the  outputs 
of  interest,  given  a  set  of  inputs. 


5.5.5  Verification  with  3D  CFD 

In  order  to  gain  confidence  in  the  proposed  approach,  it  is  imperative  that  we 
verify  the  results  obtained  using  the  Hybrid  approach  with  those  obtained  using 
detailed  3D  CFD  analysis.  A  suitably  long  segment  of  the  chosen  heat  exchanger 
geometry  was  analyzed  using  3D  CFD  method  for  a  range  of  heat  exchanger 
dimensions.  The  length  of  the  HX  segment  is  chosen  such  that  it  is  greater  than  the 
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length  of  the  entrance  region.  The  length  of  the  entrance  region  can  be  calculated 
based  on  analytical  expressions  available  in  the  literature.  The  entrance  region  needs 
to  be  accounted  for  carefully,  since  the  heat  transfer  and  pressure  drop  in  the  entrance 
region  is  very  different  than  the  one  in  the  fully  developed  flow  region.  The  same 
segments  (i.e.  same  length  and  flow  conditions)  were  analyzed  using  the  new  hybrid 
approach. 

The  lower  and  upper  bounds  for  the  different  design  variables  were  fixed  as 
shown  in  Table  5.1.  Based  on  these  bounds,  a  test  matrix  was  generated  using  a  2- 
level  full  factorial  center  augmented  design  (Myers  and  Montgomery,  2002).  Each  of 
these  designs  is  referred  to  as  a  case  for  identification  purposes.  The  3D  CFD  as  well 
as  the  hybrid  approach  was  applied  to  these  heat  exchanger  elements.  The  inlet  flow 
conditions  (velocity,  pressure  and  temperature)  were  fixed.  The  results  are  compared 
based  on  the  predicted  air  and  water  side  pressure  drop  and  the  difference  in  the  air 
and  water  temperature  across  the  heat  exchanger  are  shown  in  Figures  5.8,  5.9  and 
5.10. 

Figure  5.8  shows  the  error  (difference  between  3D  CFD  and  the  Hybrid  solver)  in 
the  predicted  air  and  water  side  pressure  drops.  The  maximum  absolute  percent  error 
for  air  side  pressure  drop  is  23%  and  that  for  water  is  45%. 
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Figure  5.8  Air  and  water  pressure  drop  errors 


Figure  5.9  shows  the  error  in  predicted  outlet  temperatures  for  air  and  water.  The 
maximum  absolute  error  for  air  outlet  temperature  is  3.5K  while  that  for  water  is 
1.6K.  Note  that  the  water  side  heat  transfer  coefficient  and  the  pressure  drop  are 
calculated  using  the  correlations  in  CoilDesigner,  namely,  the  Gnielinksi  correlation 
(Gnielinski,  1976)  for  heat  transfer  and  Churchill  correlation  (Churchill,  1977)  for 
pressure  drop. 
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Figure  5.9  Error  in  predicted  outlet  temperature 


Figure  5.10  shows  the  same  outlet  temperature  errors,  but  plotted  against  the 
water  side  Reynolds  number  (Re).  We  observe  that  more  of  the  high  errors  occur  for 
cases  with  Re  <  100.  This  brings  to  light  the  fact  that  the  correlations  in  CoilDesigner 
are  not  applicable  for  flows  with  low  Reynolds  number,  pointing  to  the  need  for  new 
correlations  for  use  in  CoilDesigner.  The  correlations  also  do  not  account  for  the 
effect  of  the  entrance  region  since  they  were  created  for  fully  developed  flow.  This 
leads  to  the  conclusion  that  the  hybrid  approach  can  replace  3D  CFD  calculations  for 
single  phase  flows  for  Re  >  100.  This  approach  was  used  in  subsequent  optimization 
studies.  It  was  also  found  that  the  average  Re  for  heat  exchanger  designs  encountered 
during  optimization  was  greater  than  100. 
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Figure  5.10  Errors  in  outlet  temperature  vs.  Re 


5.5.6  Design  Optimization  Approach 

In  general  design  optimization  is  carried  out  via  the  coupling  of  an  optimizer  with 
a  simulation.  In  this  case  the  simulation  is  the  Hybrid  solver  proposed  in  Section 
5.5.2.  In  conventional  optimization,  we  would  couple  the  Hybrid  solver  along  with  an 
optimization  algorithm.  Based  on  the  design  variables  of  interest,  a  MOGA  is  used 
for  optimization.  The  use  of  the  Hybrid  solver  with  a  MOGA  still  makes  the  task  of 
optimization  computationally  prohibitive,  since  each  function  call  will  require  a  2D 
CFD  evaluation,  which  itself  takes  about  10-60  minutes  each.  The  MOGA  will 
invoke  this  simulation  several  thousand  times.  To  reduce  this  computational  burden,  a 
previously  developed  approximation  assisted  optimization  approach  is  used. 
Metamodels  are  developed  for  the  air  side  performance  (i.e.  heat  transfer  coefficient 
and  pressure  drop)  as  a  function  of  different  geometry  parameters  such  as  An,  As,  Vs, 
v  and  w.  The  adaptive  DOE  approach  developed  in  Chapter  3  are  used  to  build  the 
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metamodels.  These  metamodels  are  then  coupled  with  CoilDesigner  to  evaluate 
performance  of  individual  heat  exchangers. 

The  design  optimization  approach  is  shown  in  the  form  of  a  flow  chart  in  Figure 
5.1 1.  The  different  steps  are  as  follows: 

Step-1:  Start  with  concept  geometry.  In  this  chapter,  the  geometry  shown  in 
Figure  5.4  is  used  as  the  starting  geometry. 

Step-2:  Parameterize  Geometry.  The  important  parameters  that  affect  the  air  side 
performance  of  the  geometry  are  identified.  This  is  done  by  conducting  a  full- 
factorial  analysis  for  all  parameters  of  interest  followed  by  analysis  using  scatter 
matrix  plots. 

Step-3:  Approximation  using  PPCFD.  In  this  step,  metamodels  are  developed  for 
air  side  heat  transfer  coefficient  and  pressure  drop  using  the  adaptive  DOE  approach 
proposed  in  Chapter  3.  The  automated  tool  PPCFD  is  used  for  conducting  the  CFD 
runs  and  post  processing.  Different  Kriging  correlations  are  evaluated  for  a  best  fit 
and  Kriging  metamodels  are  developed  for  air  HTC  and  DP. 

Step-4:  Optimization.  In  this  step,  a  MOGA  is  used  to  optimize  the  design  of  the 
heat  exchanger.  The  objective  is  to  find  optimal  heat  exchangers  which  provide  a 
given  heat  load  for  minimum  air-side  pressure  drop  and  minimum  HX  volume. 
CoilDesigner  is  used  to  assemble  heat  exchangers  and  the  metamodels  developed  in 
Step-3  are  used  to  assess  the  air  side  performance. 
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Step-1 


Step-2 


Step-3 


Figure  5.11  Proposed  design  optimization  approach  for  novel  heat  exchangers 


The  result  of  Step-4  is  a  set  of  Pareto  optimal  heat  exchanger  designs  that  deliver 
a  specified  heat  load.  The  optimization  objectives  chosen  for  analysis  in  this 
dissertation  are  air  side  pressure  drop  and  HX  volume,  but  any  other  heat  exchanger 
performance  parameter  can  be  used. 
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5.6  Application  of  Single  Response  Adaptive  DOE 


In  this  section,  the  procedure  and  the  results  for  the  application  of  single  response 
adaptive  DOE  method  (see  Chapter  3)  to  heat  exchanger  design  optimization  is 
presented. 

5.6.1  Design  Variables 

The  input  variables  used  to  build  the  metamodel  and  their  limits  are  given  in 
Table  5.1.  The  outer  diameter,  the  tube  thickness  and  the  leading  and  the  trailing 
edges  are  a  function  of  the  inner  diameter.  This  was  based  on  prior  analysis,  but  can 
be  readily  changed.  Thus  accounting  for  inner  diameter  also  accounts  for  outer 
diameter,  thickness  and  edges.  The  vertical  and  horizontal  spacing  needs  to  be 
accounted  for,  since  it  has  a  direct  influence  on  the  air-side  heat  transfer  and  pressure 
drop.  The  horizontal  spacing  also  has  an  influence  on  the  number  of  parallel  ports  a 
particular  design  will  have.  The  depth  was  chosen  such  that  the  minimum  number  of 
ports  is  1  and  maximum  number  of  ports  is  41.  Since  the  limits  imposed  on  the  inner 
diameter  differ  by  an  order  of  magnitude,  it  is  imperative  to  have  the  limits  on  the 
other  design  variables  scale  accordingly.  The  velocity  limit  was  chosen  based  on  the 
velocity  limits  for  conventional  air-conditioning  applications.  All  design  variables 
mentioned  below  are  normalized  to  the  interval  [0,  1]  when  used  in  DOE  and 
metamodel  development. 
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Table  5.1  Design  variables  for  metamodel 


Design  Variable 

Lower  Limit 

Upper  Limit 

Inner  diameter,  An 

0.1mm 

1mm 

Horizontal  spacing,  Hs  [mm] 

1.5  xfoul 

6.0  x  Aut 

Vertical  spacing,  Vs  [mm] 

2  x  Aut 

4  x  Am 

Depth,  w  [mm] 

Function  of  An 

Function  of  An 

Air  velocity,  v  [m/s] 

1.0 

3.0 

5.6.2  Metamodel  Development  -  SFCVT  Method 

The  SFCVT  method  with  the  SO  formulation  presented  in  Eq.  (3.9)  was  used.  The 

DOE  was  carried  out  separately  for  the  two  responses  viz.,  Air  DP  and  Air  HTC. 
Based  on  preliminary  approximation  studies,  it  was  found  that  approximately  200 
points  were  needed  to  achieve  MAS  of  more  than  75  with  an  acceptability  threshold 
of  10%.  The  initial  design  comprising  of  100  points  was  generated  using  the  MED 
method  and  the  SFCVT  method  was  used  to  sample  100  additional  points.  Kriging 
with  logarithmic  response  (to  avoid  negative  values  during  prediction)  was  used  to 
develop  the  metamodel  and  a  set  of  250  randomly  generated  test  points  was  used  to 
verify  the  metamodel.  The  results  of  this  verification  are  shown  in  Figures  5.12  and 
5.13. 

Figure  5.12a  shows  the  RMSE,  MAE  and  MAS  for  Air  DP  and  HTC  and  Figure 
5.12b  shows  the  RRMSE  and  RMAE  for  Air  DP  and  HTC.  It  can  be  observed  that 
more  than  80%  of  the  test  points  were  predicted  within  10%  of  the  true  value.  From 
Figure  5. 12b,  we  see  that  the  maximum  absolute  prediction  error  was  23%  for  Air  DP 
and  33%  for  Air  HTC.  These  large  errors  were  for  cases  for  which  the  CFD 
simulation  did  not  converge. 

Figure  5.13,  shows  the  box  plot  for  the  verification  errors,  and  we  see  that  the 
median  errors  values  are  close  to  0.  Based  on  these  verification  results,  the  resulting 
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metamodels  were  found  to  be  acceptable.  These  metamodels  were  then  used  for 
optimization. 


(a)  (b) 

Figure  5.12  Metamodel  verification  results  for  Air  DP  and  HTC 


Air  DP 


Air  HTC 


Figure  5.13  Metamodel  verification  for  Air  DP  and  HTC,  relative  errors 
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5.6.3  Metamodel  Development-  MSFCVT  Method 


Multiresponse  metamodels  were  also  developed  using  the  MSFCVT  method 
presented  in  Chapter  4.  The  MSFCVT  method  with  the  formulation  presented  in  Eq. 
(4.2)  was  used.  The  method  parameters  were  the  same  as  the  one  used  in  the 
development  of  the  SFCVT  metamodel,  i.e.,  100  points  in  initial  design,  100 
adaptively  sampled  points,  with  5  new  points  selected  per  stage.  The  same  random 
sample  of  250  point  was  used  for  metamodel  verification.  The  results  are  shown  in 
Figure  5.14. 


(a)  (b) 

Figure  5.14  MSFCVT  method  verification  for  air  DP  and  HTC,  (a)  RMSE,  MAE  and 

MAS,  (b)  RRMSE  and  RMAE 

From  Figure  5.14a,  we  can  see  that  the  approximately  80%  of  the  verification 
points  were  predicted  within  ±10%  of  the  actual  value  for  both  air  FITC  and  DP.  The 
maximum  absolute  error  in  predicted  air  DP  is  27%  and  for  HTC  is  36%.  These  errors 
occurred  for  cases  that  had  convergence  issues  during  CFD  simulations.  Comparing 
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the  metamodel  performance  with  that  obtained  from  SFCVT  method,  we  observe  that 
the  MSFCVT  method  yields  similarly  performing  metamodels  but  with  a  fewer 
number  of  total  function  calls,  i.e.,  300  vs.  200. 

5.7  Design  Optimization 


5.7.1  Problem  Formulation 

Based  on  the  performance  evaluation  criteria  discussed  earlier,  the  objective  of 
the  design  problem  was  to  find  new  heat  exchanger  design  which  has  minimum  air- 
side  pressure  drop  and  minimum  volume,  but  still  has  the  same  heat  duty/capacity  as 
a  1  kW  air-to-water  micro  channel  heat  exchanger. 

For  a  fair  comparison,  the  refrigerant  side  pressure  drop  (RDP)  was  constrained  to 
be  1.5  times  that  of  the  microchannel  heat  exchanger  used  in  the  comparison.  The 
optimization  problem  can  be  mathematically  represented  as: 


min  ADP,  V 

HX 

s.t. 

0  =  lkW 

ADP  <  ADPmax 

RDP  <  RDPmax 

HX  =  [Nt,  Din ,  Hs,  Vs,  w,  v] 


(5.1) 


where,  ADP  is  the  air-side  pressure  drop,  V  is  the  HX  volume,  Q  is  the  heat  load, 
ADP  and  RDP  are  the  air  and  the  refrigerant  side  pressure  drops  (respectively), 
ADPmax  and  RDPmax  are  the  limits  on  ADP  and  RDP.  Note  that  Nt  is  the  number  of 
parallel  tubes  in  the  heat  exchanger  and  is  discrete  variable,  while  An,  Hs,  Vs,  w,  v  are 
all  continuous.  The  MOOP  represented  by  Equation  5.1  was  solved  using  a  MOGA. 
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The  relevant  problem  parameters  and  constraint  limits  are  shown  in  Table  5.2.  The 
parameters  used  for  MOGA  are  given  in  Table  5.3. 


Table  5.2  Optimization  problem  parameters 


Parameter 

Value 

Heat  Load  Q,  minimum 

lkW 

Heat  Load  Q,  maximum 

1.1  kW 

ADPmax 

100  Pa 

RDPmax 

lkPa 

Air  inlet  temperature 

300K 

Water  inlet  temperature 

350K 

Air  flow  rate 

0.025  kg/s 

Water  flow  rate 

0.025  kg/s 

Table  5.3  MOGA  Parameters 


Parameter 

Value 

GA  Population  Size 

100 

GA  Population  Replacement 

10 

GA  Max  Iterations 

500 

Continuous  Variable  bit  length 

10 

The  metamodels  obtained  from  the  SFCVT  method  for  air  HTC  and  DP, 
described  in  Section  5.6.2,  were  used  in  the  CoilDesigner  calculations  for  predicting 
the  coil  performance.  The  resulting  Pareto  set  comprising  of  60  designs  is  shown  in 
Figure  5.15.  It  can  be  observed  that  all  the  Pareto  solutions  have  an  air  side  pressure 
drop  less  than  that  of  the  baseline  microchannel  coil.  There  are  37  designs  better  than 
the  baseline  microchannel  in  terms  of  HX  volume  with  the  best  heat  exchanger 
having  a  44%  less  volume  than  the  baseline  HX.  All  these  heat  exchangers  have  a 
heat  load  of  1  kW  within  +10%.  The  better  designs  also  have  an  average  of  61%  less 
material  volume  and  hence  weight. 
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Figure  5.15  Pareto  solutions  for  NGHX  approximation  assisted  optimization 

The  details  of  three  sample  designs  are  shown  in  Table  5.4,  along  with  HX 
volume,  air  pressure  drop  and  material  volume.  The  designs  shown  are  the  ones 
which  are  better  than  the  baseline  microchannel  coil.  We  observe  that  the  new  designs 
offer  a  significant  reduction  in  heat  exchanger  volume  and  the  air  side  pressure  drop 
compared  to  baseline  coil  for  the  same  application  performance. 


Table  5.4  Sample  Pareto  solutions 


Pareto 

Design 

Normalized  Design  Variables 

HX 

Vol. 

[cm3] 

Air 

DP 

[Pa] 

Matl. 

Vol. 

[cm3] 

Nt 

An 

As 

Vs 

w 

V 

1 

162 

0.316 

0.265 

0.0039 

0.180 

0.749 

126.6 

44.2 

14.52 

2 

176 

0.344 

0.252 

0.0127 

0.164 

0.427 

156.1 

23.5 

17.79 

3 

188 

0.367 

0.020 

0.0918 

0.164 

0.003 

201.6 

6.6 

22.72 

Baseline 

Micro¬ 

Channel 

NA 

230.0 

63.8 

37.3 
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Figure  5.16  shows  the  data  represented  in  Table  5.4.  The  three  sample  Pareto 
designs  are  shown  in  the  main  plot  along  with  their  relative  volume  (inset  figure) 
compared  to  the  baseline  microchannel  coil  volume. 


0  100  200  300  400  500  600  700 

1 1 X  Volume  [cm3] 


Figure  5.16  Sample  Pareto  designs  and  their  relative  volume  compared  to  the  baseline 

microchannel  coil  volume 


5.7.2  Pareto  Optimal  Designs  -  Approximation  Perspective 

In  this  section  we  analyze  the  Pareto  optimal  results  to  evaluate  how  good  the 

results  obtained  using  approximation  assisted  optimizations  are. 

The  above  optimization  was  carried  out  using  metamodels  for  the  air  side 
performance  of  the  HX  segment.  It  is  imperative  that  for  all  Pareto  solutions  we 
verify  the  predicted  air-side  performance  with  actual  CFD  runs  and  then  re-calculate 
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the  overall  performance  of  the  heat  exchanger.  The  resulting  percentage  errors  in  the 
predicted  air  side  pressure  drop  and  heat  transfer  coefficients  and  subsequently  the 
calculated  heat  load  are  shown  in  Figure  5.17.  Based  on  the  median  values,  the  air 
side  pressure  drop  is  over  predicted  by  2.8%,  the  heat  transfer  coefficient  is  under¬ 
predicted  by  1.5%  and  subsequently  the  heat  load  is  under-predicted  by  0.91%.  From 
a  heat  exchanger  perspective,  the  over-prediction  of  pressure  drop  and  under¬ 
prediction  of  HTC  is  favorable,  since  it  will  give  conservative  designs.  As  seen  from 
this  figure,  there  are  some  outliers  for  heat  load  error  values,  from  which  6  designs 
violate  the  heat  load  constraint  (i.e.  1000W  <  Q  <  1100W).  5  out  of  these  6  points 
have  a  heat  load  of  less  than  lkW  with  maximum  deviation  of  2.3%  and  one  point  has 
heat  load  more  than  1100W  with  maximum  error  of  0.5%  which  is  practically 
insignificant.  Nevertheless,  we  update  the  Pareto  set  with  the  true  responses  of  Air 
HTC  and  DP.  The  two  Pareto  sets  are  shown  in  Figure  5.18  below.  They  almost 
overlap. 
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Figure  5.17  CFD  vs.  Metamodel,  percentage  errors  in  predicted  performance 


Figure  5.18  Pareto  Designs  obtained  using  approximation  and  verification 
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5.7.3  Pareto  Optimal  Designs  -  Thermodynamic  Perspective 

In  this  section  we  provide  insight  into  the  Pareto  optimal  design  from  a 

thermodynamic  perspective. 

Careful  analysis  of  the  design  variables  corresponding  to  the  Pareto  solutions 
reveals  that  all  of  the  solutions  have  the  same  or  close  enough  normalized  inner 
diameters  in  the  range  [0.3,0.40].  This  is  because  the  inner  diameter  has  a  direct 
impact  on  the  pressure  drop  of  the  refrigerant.  MOGA  is  forced  to  drive  towards  the 
lowest  diameter  which  still  satisfies  the  pressure  drop  constraint,  while  other  variables 
such  as  the  horizontal  and  vertical  spacing  and  velocity  take  a  wide  range  of  values. 
The  HX  normalized  depth  for  the  Pareto  designs  was  also  in  a  narrow  range  of  [0.16, 
0.31],  This  is  again  due  to  the  air  side  pressure  drop  constraint.  The  designs 
corresponding  to  smallest  vertical  spacing  also  had  the  corresponding  smaller  values 
for  velocity,  which  again  is  influenced  by  the  air  side  pressure  drop  objective. 

Overall,  we  observe  that  the  new  heat  exchanger  geometry  is  capable  of 
performing  better  than  the  current  state  of  the  art  microchannel  (i.e.,  baseline 
microchannel  coil).  This  improvement  is  mainly  due  to  significant  compactness  of  the 
new  geometry.  In  other  words,  the  new  geometry  can  accommodate  more  heat 
transfer  area  per  unit  HX  volume  and  per  unit  HX  material  volume.  The  above 
comparison  was  carried  out  for  one  application,  but  using  the  newly  developed 
approach  optimal  HX  designs  for  other  applications  (i.e.,  different  fluid  inlet 
conditions  and  heat  load)  can  be  readily  obtained.  As  a  final  note,  it  should  be  pointed 
out  that  the  above  results  should  be  validated  experimentally. 
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5.8  Summary 


In  this  chapter,  a  new  approach  for  the  multi-level  performance  evaluation  and 
design  optimization  of  novel  heat  exchanger  elements  is  introduced.  The  objective  of 
developing  this  approach  is  to  allow  streamlined  search  of  novel  heat  exchanger 
geometries  for  air-to-refrigerant  heat  exchanger  applications.  Instead  of  conducting 
full  fledged  3D  CFD  simulations,  the  new  approach  uses  2D  CFD  techniques  coupled 
with  an  s-NTU  solver  for  heat  exchanger  performance  evaluation.  This  results  in  over 
95%  savings  in  computation  time.  The  2D  CFD  is  used  for  air-side  performance 
evaluation  of  small  heat  exchanger  segment.  The  s-NTU  solver  is  then  used  to 
assemble  these  small  HX  segments  or  elements  to  design  full  scale  heat  exchangers. 
The  new  approach  was  verified  for  a  range  of  dimensions  for  a  heat  exchanger.  The 
comparison  of  the  results  (viz.,  pressure  drop  and  change  in  temperature)  showed  that 
the  new  hybrid  approach  is  valid  for  applications  with  single  phase  flows  with 
Reynolds  number  greater  than  100. 

Even  though,  the  replacement  of  3D  CFD  with  the  new  approach  results  in  95% 
savings  of  computational  effort,  the  execution  time  is  still  prohibitive  (~  1  week  for 
one  optimization  run).  We  apply  the  single  response  adaptive  DOE  approach 
developed  in  Chapter  3  to  this  problem,  where  in  we  approximate  the  air  performance 
(DP  and  HTC)  obtained  from  the  2D  CFD  calculations.  The  resulting  metamodels  for 
air  side  DP  and  HTC  had  MAS  of  80  which  was  acceptable  for  the  current 
optimization  problem.  The  fact  that  MAS  of  80  is  acceptable  is  further  reinforced  by 
the  post-optimization  verification  results. 


168 


The  optimization  problem  was  formulated  as  an  MOOP,  minimizing  the  heat 
exchanger  volume  and  the  air  side  pressure  drop  while  providing  a  heat  load  of  1  kW. 
A  conventional  state-of-the  art  microchannel  heat  exchanger  was  used  as  a  baseline 
for  comparison.  The  Pareto  optimal  designs  resulting  from  the  MOOP  solution  had  37 
heat  exchanger  designs  that  are  better  than  the  baseline  heat  exchanger,  with  a  volume 
reduction  of  44%  and  material  reduction  of  61%.  In  addition,  the  computational 
savings  are  enormous,  since  the  metamodels  can  now  be  used  for  additional 
optimization  studies  with  varying  constraints  etc.  The  Pareto  designs  obtained  from 
the  approximation  assisted  optimization  were  also  verified  against  full  CFD 
calculations.  It  was  found  that  amongst  the  Pareto  solutions  obtained  using 
approximation,  6  designs  were  infeasible  but  by  a  very  marginal  amount,  while  all 
other  designs  overlap  with  those  obtained  after  CFD  verification. 

The  inner  diameter  and  the  depth  of  the  heat  exchangers  corresponding  to  the 
Pareto  solutions  were  distributed  in  a  very  narrow  range.  The  refrigerant  side  pressure 
drop  constraint  was  the  dominating  factor  in  the  choice  of  inner  diameters,  while  the 
air  side  pressure  drop  objective  was  the  dominating  factor  for  the  depth  of  the  heat 
exchanger.  This  is  because,  small  inner  diameter  leads  to  higher  refrigerant  side 
pressure  drop  and  larger  depth  leads  to  larger  air  side  pressure  drop. 
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Chapter  6  Conclusion 


In  this  dissertation  three  research  thrusts  are  presented.  The  research  thrusts  are, 

(a)  adaptive  DOE  approach  for  single  response  deterministic  computer  simulations, 

(b)  adaptive  DOE  approach  for  multiresponse  deterministic  computer  simulations  and 

(c)  approximation  based  framework  for  the  performance  evaluation  and  optimization 
of  air-cooled  heat  exchangers. 

This  chapter  is  organized  as  follows:  Section  6.1  provides  a  summary  of  the  three 
research  thrusts  followed  by  conclusions  in  Section  6.2.  Section  6.3  highlights  the 
main  contributions  of  this  dissertation  and  section  6.4  provides  a  discussion  for 
potential  future  research  direction. 

6.1  Summary 

In  this  section,  the  three  research  thrusts  are  summarized  as  follows: 

(a)  Adaptive  DOE  for  Single  Response:  In  Chapter  3  an  adaptive  DOE  method 
for  single  response  detenninistic  computer  experiments  is  presented.  The 
proposed  approach  is  based  on  the  notion  of  tradeoff  that  needs  to  be  made 
by  a  DOE  method  between  allocating  samples  in  sparsely  populated 
regions  of  the  design  space  versus  allocating  samples  in  regions  in  which 
the  response  is  most  sensitive.  The  sensitivity  of  the  response  in  a 
particular  region  of  the  design  space  is  gauged  by  the  calculation  of  leave 
one  out  cross-validation  error,  while  the  sparsely  sampled  regions  are 
tracked  through  the  use  of  a  space  -filling  criterion.  The  approach  is  such 
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that  one  point  is  sampled  per  stage,  though  there  are  provisions  for 
sampling  multiple  points  per  stage  to  take  advantage  of  parallel  computing 
capabilities.  The  proposed  approach  is  then  applied  to  a  suite  of  test 
problems  compiled  as  a  part  of  this  dissertation. 

(b)  Adaptive  DOE  for  Multiple  Responses:  Chapter  4  extends  the  approach 
proposed  in  Chapter  3  to  multiple  responses.  The  problem  of  adaptive 
DOE  for  multiresponse  case  is  posed  as  a  multiobjective  optimization 
problem  where  the  tradeoff  is  between  allocating  new  samples  such  that 
they  simultaneously  improve  the  metamodel  performance  for  each  of  the 
responses.  A  multiobjective  genetic  algorithm  is  used  to  solve  the 
optimization  problem  resulting  in  a  set  of  Pareto  optimal  solutions.  A  level 
of  dominance  criteria  is  introduced  to  select  one  point  from  the  Pareto  set. 

(c)  Approximation  Based  Framework  for  Design  of  Novel  Heat  Exchangers: 
In  Chapter  5  a  new  approximation  based  framework  is  proposed  to 
evaluate  the  multi-level  performance  of  air-cooled  heat  exchangers.  An 
approach  that  uses  2D  CFD  coupled  with  an  s-NTU  heat  exchanger  model 
is  developed  and  tenned  as  the  Hybrid  solver.  This  Hybrid  solver  allows  us 
to  achieve  the  same  results  as  those  obtained  using  a  detailed  3D  CFD 
analysis  but  with  an  order  of  magnitude  reduction  in  computational  cost. 
The  performance  evaluation  approach  is  then  coupled  with  a  multi¬ 
objective  genetic  algorithm  to  find  optimal  designs  for  novel  air-cooled 
heat  exchanger. 


171 


6.2  Conclusions 


This  section  provides  some  concluding  remarks  on  each  research  thrust. 

6.2.1  Single  Response  Adaptive  DOE 

The  cross  validation  based  single  response  adaptive  DOE  approach  outperforms 
existing  adaptive  methods  such  as  MSD  (Jin,  et  ah,  2002),  ACE  (Li  and  Azarm,  2006; 
Li,  2007)  and  non-adaptive  method  such  as  the  MED  (Shewry  and  Wynn,  1987; 
Koehler  and  Owen,  1996)  method  in  more  than  50%  of  the  test  problems.  The 
comparison  was  based  on  the  median  values  of  the  metamodel  performance  metrics. 
This  observation  is  supported  by  the  descriptive  statistics  drawn  from  the  results 
generated  by  applying  the  said  approaches  to  the  24  test  problems.  In  certain  cases  it 
was  observed  that  even  though  the  proposed  approach  does  not  perfonn  better  in 
terms  of  RMSE  and  MAE  in  comparison  to  other  approaches,  it  does  perform  well  in 
terms  of  the  relative  predictions,  i.e.,  RRMSE  and  RMAE.  Even  though  RMSE  and 
MAE  are  widely  used  in  the  literature  for  comparing  metamodel  performance,  the 
results  in  this  dissertation  show  that  these  measures  are  not  enough  to  judge  a 
metamodel  and  as  such  additional  metric  such  MAS  is  also  used.  The  metrics  of 
RMSE,  MAE,  RRMSE,  RMAE  and  MAS  are  used  in  this  dissertation  to  compare  the 
performance  of  different  metamodels. 

As  detailed  in  Chapter  3,  the  proposed  approach  develops  a  metamodel  of  LOO 
errors  to  predict  the  regions  of  high  sensitivity  in  the  response  space.  For  certain 
problems  (e.g.,  test  problem  P10  in  Chapter  3),  these  intermediate  metamodels  may 
not  be  accurate  enough  thus  resulting  in  poor  sampling  which  in  turn  leads  to 
underperforming  (compared  to  non-adaptive  method)  metamodels.  It  is  important  to 
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objectively  evaluate  the  intermediate  metamodels  at  each  stage  and  if  the  metamodel 
is  not  acceptable,  then  the  adaptive  method  must  switch  to  a  purely  space-filling 
method  such  as  MED.  It  should  be  pointed  out  that  the  proposed  approach  aims  at 
sampling  for  a  globally  accurate  metamodel.  Thus,  it  yields  better  metamodels  for 
problems  with  non-symmetric  CAMM  regions,  where  as  non-adaptive  space-filling 
methods  are  better  suited  for  test  problems  with  symmetric  and  highly  distributed 
CAMM  regions. 

In  all  the  numerical  experiments,  only  the  Gaussian  correlation  was  used  for  the 
Kriging  metamodels.  For  certain  types  of  problems,  this  correlation  may  not  be  the 
best  correlation  and  hence  during  implementation  multiple  correlations  need  to  be 
evaluated  and  the  best  chosen. 

6.2.2  Multiresponse  Adaptive  DOE 

The  cross  validation  based  multiresponse  adaptive  DOE  approach  outperforms 
existing  adaptive  DOE  techniques  (i.e.,  DMM  method  by  Li,  et  ah,  2006,  MCVV 
method  by  Romero,  et  ah,  2006)  based  on  the  results  published  in  the  literature.  The 
newly  developed  approach  was  also  applied  to  a  suite  of  8  test  problems  with  multiple 
inputs  and  multiple  responses.  Based  on  the  numerical  experiments,  it  was  found  that 
the  MSFCVT  approaches  yields  comparable  or  better  metamodels  for  all  responses 
simultaneously  when  compared  with  those  obtained  using  the  MED  (Shewry  and 
Wynn,  1987;  Koehler  and  Owen,  1996)  method. 

The  conclusions  drawn  for  the  single  response  adaptive  DOE  method  in  Section 
6.2.1  are  also  applicable  to  the  multiresponse  adaptive  DOE  method. 
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It  was  observed  at  each  stage,  that  MSFCVT  simultaneously  tried  to  improve  the 
metamodels  for  each  of  the  multiple  responses.  This  can  sometimes  lead  to  over¬ 
improvement  in  some  of  the  metamodels  as  observed  for  Numerical  Example- 1.  At 
some  point  during  the  adaptive  DOE  phase,  the  performance  of  individual 
metamodels  needs  to  be  evaluated  and  the  MSFCVT  method  needs  to  focus  on  the 
responses  with  less  accurate  metamodels.  This  cannot  be  carried  out  in  the  current 
implementation,  since  no  assumption  about  availability  of  a  test  sample  for 
verification  was  made. 

6.2.3  Approximation  Based  Framework  for  Design  of  Novel  Heat  Exchangers 

A  new  approach  for  multi-level  performance  evaluation  of  novel  heat  exchanger 

geometries  that  uses  approximation  assisted  optimization  is  developed  and  applied  to 
a  real  world  case  study. 

The  proposed  Hybrid  approach  is  valid  for  single  phase  flows  with  Re  >  100. 
This  is  because  the  correlations  in  CoilDesigner  are  not  valid  in  the  low  Re  range. 
This  results  in  large  calculation  errors  in  the  water  side  heat  transfer  coefficients  and 
pressure  drop,  which  in  turn  affects  the  air-side  performance,  leading  to  large 
discrepancies  between  the  heat  exchanger  performance  obtained  using  a  full  3D  CFD 
simulation  and  those  obtained  using  the  hybrid  approach. 

In  the  study,  individual  metamodels  were  developed  for  the  air  side  DP  and  HTC, 
with  100  initial  points  and  100  new  points  for  each  metamodel,  resulting  in  a  total  of 
300  function  evaluations  (i.e.,  initial  100  points  were  common).  The  metamodels 
were  verified  using  250  randomly  generated  samples.  The  resulting  MAS  was  80.0 
for  both  air  HTC  and  DP.  Optimization  was  carried  out  with  these  metamodels  and 
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the  results  were  verified  against  actual  CFD  runs.  The  comparison  was  surprisingly 
good  with  the  median  prediction  errors  of  2.8%  (over-prediction),  1.5%  (under- 
prediction)  and  0.9%  (under-prediction)  for  air  DP,  HTC  and  heat  load  respectively, 
leading  to  the  conclusion  that  for  multi-level  analysis  (e.g.,  simulation  at  heat 
exchanger  segment  level  and  simulation  at  the  heat  exchanger  level)  a  low  fidelity, 
i.e.,  relatively  less  accurate  metamodel  may  be  acceptable.  This  is  because  at  the  top 
level,  the  errors  in  low  level  metamodels  do  not  have  a  significant  impact.  If  an 
uncertainty  analysis  can  be  carried  out  in  early  stages  of  analysis  and  optimization,  it 
could  result  in  computational  savings,  since  the  metamodels  need  not  be  highly 
accurate. 

Analysis  of  the  Pareto  optimal  designs  led  to  the  conclusion  that  the  optimization 
search  was  dominated  by  the  refrigerant  side  pressure  drop  constraint  which  is  a 
function  of  the  tube  inner  diameter.  All  the  Pareto  solutions  had  nonnalized  Dm  in  the 
interval  [0.3,  0.4],  Another  factor  heavily  driving  the  search  was  the  air-side  pressure 
drop  objective,  which  is  influenced  by  the  depth  of  the  heat  exchanger  in  the  air-flow 
direction  in  addition  to  air  velocity  and  vertical  spacing.  All  Pareto  designs  had 
normalized  depth  in  the  range  [0. 16,  0.3 1]. 

Based  on  the  promising  designs  obtained  using  the  proposed  approach,  it  is 
concluded  that  this  tool  will  significantly  improve  the  way  novel  heat  exchangers 
with  the  proposed  and  similar  geometry  are  explored. 
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6.3  Contributions 


The  main  research  contributions  of  this  dissertation  are: 

1.  An  approach  for  single  response  adaptive  DOE  for  deterministic  computer 
simulations  has  been  developed.  A  new  problem  formulation  is  proposed  that 
views  adaptive  DOE  problem  as  a  multiobjective  optimization  problem. 
Applicability  of  the  proposed  approach  has  been  tested  on  a  broad  set  of  test 
problems  from  the  literature.  It  is  shown  that  on  an  average  the  proposed 
approach  yields  better  metamodel  for  the  same  number  of  function  calls,  than 
the  approaches  in  the  literature. 

2.  An  approach  for  multiresponse  adaptive  DOE  has  been  developed.  A  new 
multiobjective  problem  formulation  is  proposed  which  makes  tradeoffs 
between  the  multiple  responses.  Applicability  of  the  proposed  approach  is 
verified  by  successfully  applying  it  to  the  test  problems  from  literature.  Based 
on  comparison  with  published  results,  it  is  shown  that  the  proposed  approach 
performs  better  than  two  other  adaptive  DOE  methods  from  the  literature.  It  is 
shown  that  in  multiresponse  DOE,  there  is  a  tradeoff  between  allocating 
points  to  improve  one  response  vs.  another. 

3.  A  new  multi-level  performance  evaluation  and  optimization  approach  has 
been  developed  for  design  of  novel  air-cooled  heat  exchangers.  The  approach 
helps  reduce  the  computation  time  from  a  few  months  to  a  couple  of  days, 
resulting  in  over  95%  savings  in  computational  effort  without  any  practical 
loss  of  accuracy  in  the  results.  In  this  approach,  the  heat  exchanger 
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performance  is  evaluated  not  only  at  the  element  level,  which  is  the 
conventional  approach,  but  also  at  the  component  level.  The  component  level 
evaluation  tells  us  upfront  as  to  whether  the  new  geometry  can  compete  with 
current  state  of  the  art  heat  exchangers.  The  multi-level  performance 
evaluation  approach  along  with  the  approach  developed  in  1  above  is  used  in 
the  development  of  a  multi-objective  genetic  algorithms  based  optimization 
tool  for  the  optimization  of  novel  heat  exchanger.  The  new  heat  exchanger 
design  has  44%  less  volume  and  uses  61%  less  material  compared  to  current 
state  of  the  art  microchannel  heat  exchangers.  The  new  heat  exchanger  design 
approach  is  not  limited  to  conventional  geometries  (e.g.  round  tubes)  and  can 
be  applied  to  non-conventional  heat  exchanger  geometries  such  as  tubes  with 
aero-foil  shaped  cross  section. 

6.4  Future  Research  Directions 

The  research  presented  in  this  dissertation  has  addressed  several  shortcomings  of 
adaptive  DOE  approaches  proposed  in  the  literature.  In  order  to  make  the  approach 
more  efficient  and  less  sensitive  to  initial  samples,  the  following  research  directions 
could  be  pursued: 

(a)  The  examples  shown  in  this  dissertation  used  the  proposed  approach  to 
develop  metamodels  for  the  entire  design  space  (global  metamodels).  But  the 
approach  can  also  be  used  with  a  restricted  design  space.  One  example  of  such 
an  application  is  online  approximation,  where  the  M/SFCVT  method  can  be 
used  to  intelligently  sample  new  points  in  the  regions  of  expected  optima. 
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(b)  The  resource  allocation  problem,  wherein  the  total  number  of  available 
function  calls  is  fixed  and  then  the  approach  should  be  able  to  judiciously  and 
objectively  allocate  samples  for  initial  design,  subsequent  sequential  sampling 
and  for  random  verification.  Moreover,  the  approach  needs  to  handle  this 
without  any  prior  knowledge  of  the  computer  simulation.  A  preliminary 
approach  is  proposed  in  Section  3.6.  But  a  clever  way  to  evaluating 
metamodel  performance  and  stalling  condition  without  additional  function 
calls  is  required. 

(c)  In  developing  the  metamodels  for  the  LOO  errors,  a  Gaussian  correlation  was 
used.  Other  correlation  options  need  to  be  investigated  along  with  any 
additional  tuning.  At  each  stage  the  best  LOO  model  needs  to  be  developed 
and  since  no  information  about  the  simulation  is  known,  this  points  to  a  trial 
and  error  strategy  where  each  type  of  correlation  could  be  tried  and  compared. 

(d)  In  the  studies  presented,  the  MED  method  was  used  to  generate  all  initial 
designs.  Other  space  filling  techniques  such  as  Latin-hypercubes  should  be 
investigated. 

(e)  The  proposed  cross-validation  approach  can  be  extended  to  non-Kriging  based 
metamodels. 

(f)  Develop  a  convenient  graphical  user  interface  for  the  approximation  based 
framework  for  performance  evaluation  of  heat  exchangers,  which  will  allow 
end-users  from  the  industry  to  quickly  and  efficiently  evaluate  novel  heat 
exchanger  geometries. 
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(g)  In  Chapter  5,  for  the  geometry  under  consideration,  no  significant  structural 
analysis  was  required.  But  in  real  life,  new  concepts  may  require  careful 
tradeoff  between  structural  aspect  (i.e.,  material  used)  and  the  heat  transfer 
performance  giving  rise  to  a  multidisciplinary  optimization  problem  possibly 
a  coupled  problem.  A  possible  multiobjective  collaborative  optimization 
approach  such  as  the  one  presented  in  Aute  and  Azarm  (2006)  could  be 
implemented  as  a  part  of  the  framework.  The  use  of  approximation  is  a  must 
in  such  problems  since  the  computational  effort  increases  super  linearly  with 
the  number  of  different  disciplines  involved. 
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Appendix-A  Test  Problem  Suite 


In  order  to  evaluate  the  general  applicability  of  the  approaches  proposed  in  this 
dissertation,  a  suite  of  test  problems  was  compiled.  Some  of  these  test  problems  are 
based  on  those  found  in  the  literature  and  used  by  other  researchers  (Jin  et  ah,  2002; 
Li  and  Azarm,  2006;  Romero  et  ah,  2006),  while  others  were  devised  as  a  part  of  this 
dissertation  (see  below).  The  test  problems  are  used  to  evaluate  the  following  aspects 
of  an  adaptive  design  of  experiments  approach:  (a)  scalability  to  more  than  two 
dimensions,  and  (b)  applicability  to  multi-modal  and  highly  sensitive  responses.  The 
remainder  of  this  Appendix  provides  short  background  and  the  literature  references 
for  each  of  the  test  problems,  followed  by  the  equations  and  the  domain  of  each 
problem. 

Problem  PI  and  P4  are  from  Dixon  and  Szego  (1978),  P2,  P13,  P14  and  P16  are 
made  up  test  problems.  P5  is  Griewank’s  function  (Hedar,  2005),  P6  is  Himmelblau’s 
function  and  P7  is  the  hump  function  (Floudas  et  ah,  1999),  Pll  is  Rastrigin’s 
function  (Hedar,  2005)  ,  P12  is  SchwefeTs  function  (Hedar,  2005)  and  P17  is 
Coleville  function  (Hedar,  2005),  P15  is  the  wine  bottle  function  (Van  Veldhuizen 
and  Lamont,  1998).  Problem  P9  is  Osio’s  function  (Osio  and  Amon,  1996),  P8  is 
Sasena’s  function  (Sasena  et  ah,  2000)  and  P10  is  Busby’s  (Busby  et  ah,  2007) 
function.  Problems  PI 8,  P19  and  P20  are  the  scaled  4  dimensional  versions  of  their  2 
dimensional  counterparts  while  problems  P22,  P23  and  P24  are  scaled  to  8 
dimensions,  i.e.,  8  inputs.  Problem  P21  is  the  borehole  model  simulation  from  Joseph 
et  al„  (2008) 
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Table  A.l  Single  response  test  problems 


Problem 


PI:  Bran  in  Function  (Dixon  and  Szego,  1978) 


I  x  \  5x  f  11 

y  =  [x2- 5.1  -1-  2h - ^-6]2  +  10  1 - cos(x,)  +  10 

v  2tt  )  7t  v  87Z" ) 


Xj  e[-5,10],x2  e[0, 15] 


P2:  Trigonometric  Function- 1 


y  =  cosUxf  +  x2 ' 


x1?x2  e  [-5,5] 


P3:  Dixon-Price  Function  (Hedar,  2005) 


y  =  (X  -  O2  +  Z  ~  xm  )2 

i=2 

x.  e  [-10,10],/  = 


P4:  Goldstein-Price  (Dixon  and  Szego,  1978) 

v  =  [1  +  (Xj  +x2  + 1)2  (1 9  —  14xj  -l-  3x2  -  14x2  +6xjX2  +3x2)][30  +  (2x1  -3x2)2x 
(18-32xj  +12xj2  +48x2  -6xxx2  +27x2)] 

Xj,x2  e[-2,2] 

P5:  Griewank  Function  (Hedar,  2005) 


v  =  V  — - FT cos(x,  /  Ji)  +1 

'  tf  4000  if  ‘  V 

x.  e  [-20,20],/  = 


6. 

P6:  Himmelblau’s  Function  (Floudas  et  al.,  1999) 

y  =  (xf  +  x2  - 1 1)2  +  (x,  +  x2  -  7)2 

Xj,x2  e[-4,4] 

7. 

P7:  Hump  Function  (Floudas  et  al.,  1999) 

v  =  4xj2  -  2 .  lxj4  +  ^  xf  +  XjX2  -  4x2  +  4x2 
x1?x2  e[-5,5] 

8. 

P8:  Sasena’s  Function  (Sasena  et  al.,  2000) 

v  =  2  +  0.0 l(x2  -  xf)2  +  (1  —  Xj )2  +  2(2-x2)  +  7sin(0.5x1)sin(0.7x1x2) 
x1?x2  e[0,5] 

9. 

P9:  Osio’s  Function  (Osio  and  Amon,  1996) 

v  =  cos(6(Xj  —  0.5))  +  3. 1(|  Xj  -0.7  )  +  2(xj  -0.5)  +  sin( - - - )  +  0.5x, 

Xj  -  0.5  +0.31 

XpXj  e[0,l] 

10. 

P10:  Busby’s  Function  (Busby  et  al.,  2007) 

sin(^/xj  +x2  )  +  10  1/2 

y  =  7  r- - -  +3  Xj  x2| 

V*f+*2 

Xj,x2  e  [-10,10] 

11. 

PI  1:  Rastrigin’s  Function  (Hedar,  2005) 

y  =  20  +  xj2  +  x2  - 1 0(cos(2^Xj )  +  cos(2^x, )) 

Xj,X2  €[-1,1] 
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12. 

P12:  Schwefel  Function  (Hedar,  2005) 

y  =  418.9829h  sin  J\x,  |) 

(=i 

x,.  e  [-500,500],/'  =  1 

n  =  2 

13 

P13:  Trigonometric  Function-2 

y  =  sin^  +  x2) 

Xj  e  [0,2;r],x,,  e[0,;r] 

14. 

P14:  Trigonometric  Function-3 

v  =  sin(Xj)cos(x2) 

Xj,x2  e  [-3,3] 

15. 

PI 5:  Wine  Bottle  Function  (Van  Veldhuizen  and  Lamont,  1998) 

C  +  l 
v  = 

'  9.32 

C  —  ^  +  sin(x4  +  B  +  2) 

A  =  6xj  -  3  ;  2?  =  6x2  -  3 
x1?x2  e[0,l] 

16. 

PI 6:  Numerical  Function- 1 

V  =  XjX,  /  (x2  +  x2  +  0. 1) 

Xj,x2  e  [-10,10] 
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17. 

PI 7:  Coleville  Function  (Hedar,  2005) 

a,  =  x,  -1,  i  =  1,...,4 

v  =  100(xf  —  x9 )-  +  o]"  +  <7^  +  90(x^  —  x4)  + 10.  1(<?t  +  )  + 19.8  ci~,ci^ 
x,  e  [-10,10],  i  =  1,...,4 

18. 

PI 8:  Dixon-Price  Function,  4D 

Same  as  P3  with  n=4 

19. 

PI 9.  Griewank’s  Function,  4D 

Same  as  P5,  with  n  =  4 

20. 

P20:  Schwefel  Function,  4D 

Same  as  P12  with  n=4 

21. 

P21:  Borehole  model  (Joseph  et  ah,  2008) 

2  kTu{Hu-H}) 

T  2T  T  T  1 

In  (r/r  )  1+  %  +  “ 

L  Hr/rw)rXv  7]  J 

rw  e  [0.05,0.15],  re  [100,50000] 

Tu  e  [63070,1 15600],  Hu  e  [990, 11 10] 

7]  e  [63. 1,1 16],  77,  e  [700,820] 

L  e[1120,1680],  Kw  e[9855,12045] 

22. 

P22:  Dixon-Price  Function,  8  Dimensions 

Same  as  P3,  with  n=8 

23. 

P23:  Griewank’s  Function,  8  Dimensions 

Same  as  P5,  with  n  =  8 

24. 

P24:  Scwefel’s  Function,  8D 

Same  as  P12  with  n=8 
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On  the  following  pages  the  respective  response  surfaces  are  shown  for  test 
problems  PI  through  P16  with  inputs  nonnalized  to  [0,  1],  The  normalization  is 
carried  out  by  linearly  mapping  the  ranges  of  input  variables  shown  in  Table  A.l  for 
each  input  to  the  interval  [0,  1],  As  an  example,  in  problem  P5,  the  input  range  of 
[-20,  20]  is  changed  to  [0,  1]  by  replacing  x  with  (x+20)/40.  (x+20)/40  has  the  range 
[0,  1]  as  shown  in  the  figures. 
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(a)  (b) 

Figure  A.l  Response  surface  for  problem  (a)  PI  and  (b)  P2 


(a)  (b) 

Figure  A.2  Response  surface  for  problem  (a)  P3  and  (b)  P4 
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Figure  A.3  Response  surface  for  problem  (a)  P5  and  (b)  P6 
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(a) 


(b) 


Figure  A.4  Response  surface  for  problem  (a)  P7  and  (b)  P8 
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(a) 


(b) 


Figure  A.5  Response  surface  for  problem  (a)  P9  and  (b)  P10 


(a) 


(b) 


Figure  A.6  Response  surface  for  problem  (a)  Pll  and  (b)  P12 
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(a)  (b) 

Figure  A.7  Response  surface  for  problem  (a)  P13  and  (b)  P14 


1 


(a)  (b) 

Figure  A.8  Response  surface  for  problem  (a)  P15  and  (b)  P16 
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In  Table  A. 2,  the  multiresponse  test  problems  are  described.  The  multiresponse  test 
problems  are  basically  a  combination  of  the  single  response  test  problems  used  to 
mimic  a  multiresponse  black-box  simulation.  The  test  problems  are  labeled  PM1- 
PM8.  For  each  response,  the  corresponding  problem  number  from  the  single  response 
test  suite  from  Table  A.l  is  given. 


Table  A.2  Multiresponse  test  problems 


Sr. 

No. 

Problem 

1. 

Problem  PM  1 

Vi=  P17,  V2  =  P18 

2. 

Problem  PM2 

yi=P19,y2  =  P20 

3. 

Problem  PM3 

yi=  P21,y2  =  P22 

4. 

Problem  PM4 

yi=  P23,y2  =  P24 

5. 

Problem  PM5 

y  1=  P  k  V2  =  P6,  y3=  P8,  y 4  =  P3 

6. 

Problem  PM6 

y i=  P4,  y2  =  p  1 1 ,  v3=  P7,  y4  =  P 1 0 
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